Madrid House Price
Sergio Fernández, Miguel Ocón, Enrique Roa
30 de junio de 2023
1 Introducción y definición de objetivos
El análisis de datos en el mercado inmobiliario es algo muy común y que lleva muchísimos años desarrollándose. Con el objetivo de predecir que aspectos influyen principalmente en el precio de las casas en Madrid, hemos seleccionado un dataset (kaggle - Madrid House Price) con viviendas en venta de la capital española.
library(readr)
library(dplyr)
library(tidyr)
library(kableExtra)
library(ggplot2)
library(GGally)
library(gridExtra)
library(cowplot)
library(ggcorrplot)
library(gmodels)
library(caret)
library(ggfortify)
library(scales)
library(class)
library(distances)
library(visreg)
library(rpart)
library(rpart.plot)
library(rattle)
library(randomForest)
library(e1071)
library(pROC)
library(cluster)
library(tidyverse)
library(stringr)
library(Metrics)
library(factoextra)
library(NbClust)
library(neuralnet)
# Librerías para SERIES TEMPORALES
library(xts)
library(forecast)
library(tsibble)
library(tidyquant)
library(feasts)
library(fable)
library(stats)
library(zoo)
library(DT)paste(R.Version()$version.string)## [1] "R version 4.1.2 (2021-11-01)"2 Análisis exploratorio inicial
2.1 Lectura y preparación de los datos
En esta primera etapa cargamos el dataset descargado y hacemos una observación superficial de las variables recogidas y sus características. Se representan también las 10 filas iniciales.
mhp <- read_csv("house_price_madrid_14_08_2022.csv")head(mhp, 10) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| price | house_type | house_type_2 | rooms | m2 | elevator | garage | neighborhood | district |
|---|---|---|---|---|---|---|---|---|
| 495000 | planta 1 | exterior | 3 | 118 | TRUE | TRUE | Chopera | Arganzuela |
| 485000 | planta 2 | exterior | 2 | 82 | TRUE | TRUE | Palos de Moguer | Arganzuela |
| 315000 | planta 2 | exterior | 2 | 72 | FALSE | FALSE | Legazpi | Arganzuela |
| 585000 | planta 4 | exterior | 2 | 174 | TRUE | TRUE | Palos de Moguer | Arganzuela |
| 255000 | bajo | exterior | 3 | 75 | FALSE | FALSE | Acacias | Arganzuela |
| 299000 | planta 1 | exterior | 1 | 69 | TRUE | FALSE | Chopera | Arganzuela |
| 265000 | planta 3 | exterior | 2 | 54 | TRUE | FALSE | Delicias | Arganzuela |
| 290000 | planta 3 | interior | 4 | 69 | TRUE | FALSE | Chopera | Arganzuela |
| 660220 | planta 2 | exterior | 3 | 129 | TRUE | TRUE | Imperial | Arganzuela |
| 525000 | planta 4 | exterior | 4 | 111 | TRUE | FALSE | Chopera | Arganzuela |
summary(mhp)## price house_type house_type_2 rooms
## Min. : 725 Length:15975 Length:15975 Min. : 1.000
## 1st Qu.: 195000 Class :character Class :character 1st Qu.: 2.000
## Median : 359973 Mode :character Mode :character Median : 3.000
## Mean : 624233 Mean : 2.848
## 3rd Qu.: 749000 3rd Qu.: 3.000
## Max. :13950000 Max. :41.000
## m2 elevator garage neighborhood
## Min. : 1.0 Mode :logical Mode :logical Length:15975
## 1st Qu.: 66.0 FALSE:4597 FALSE:11528 Class :character
## Median : 93.0 TRUE :11378 TRUE :4447 Mode :character
## Mean :124.8
## 3rd Qu.:142.0
## Max. :989.0
## district
## Length:15975
## Class :character
## Mode :character
##
##
## str(mhp) ## spc_tbl_ [15,975 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ price : num [1:15975] 495000 485000 315000 585000 255000 ...
## $ house_type : chr [1:15975] "planta 1" "planta 2" "planta 2" "planta 4" ...
## $ house_type_2: chr [1:15975] "exterior" "exterior" "exterior" "exterior" ...
## $ rooms : num [1:15975] 3 2 2 2 3 1 2 4 3 4 ...
## $ m2 : num [1:15975] 118 82 72 174 75 69 54 69 129 111 ...
## $ elevator : logi [1:15975] TRUE TRUE FALSE TRUE FALSE TRUE ...
## $ garage : logi [1:15975] TRUE TRUE FALSE TRUE FALSE FALSE ...
## $ neighborhood: chr [1:15975] "Chopera" "Palos de Moguer" "Legazpi" "Palos de Moguer" ...
## $ district : chr [1:15975] "Arganzuela" "Arganzuela" "Arganzuela" "Arganzuela" ...
## - attr(*, "spec")=
## .. cols(
## .. price = col_double(),
## .. house_type = col_character(),
## .. house_type_2 = col_character(),
## .. rooms = col_double(),
## .. m2 = col_double(),
## .. elevator = col_logical(),
## .. garage = col_logical(),
## .. neighborhood = col_character(),
## .. district = col_character()
## .. )
## - attr(*, "problems")=<externalptr>El dataset contiene 15.975 observaciones (correspondientes a una vivienda cada una) y 9 variables (de las cuales son 6 cualitativas y 3 cuantitativas).
A continuación, la descripción de cada una de las variables:
price: precio
house_type: tipo de vivienda (casa, chalet, piso…)
house_type_2: si es exterior o interior
rooms: número de habitaciones
m2: metros cuadrados
elevator: si tiene ascensor
garage: si incluye garaje
neighborhood: barrio de Madrid
district: distrito de Madrid
Tras esta primera observación se harán algunas modificaciones en el dataset para transformar a variables de tipo categórico las que correspondan, así como transformar a binarias (1/0) aquellas con sólo dos valores posibles.
# Preparación de los datos. Convertimos a factores las variables que son categóricas y a 1-0 las binarias.
mhp <- mhp %>% rename(exterior = house_type_2)
mhp$exterior = ifelse(mhp$exterior == "exterior", 1, 0)
mhp$exterior = factor(mhp$exterior, levels = c(1,0))
mhp$elevator = ifelse(mhp$elevator == "TRUE", 1, 0)
mhp$elevator = factor(mhp$elevator, levels = c(1,0))
mhp$garage = ifelse(mhp$garage == "TRUE", 1, 0)
mhp$garage = factor(mhp$garage, levels = c(1,0))
mhp$house_type <- as.factor(mhp$house_type)
mhp$exterior <- as.factor(mhp$exterior)
mhp$elevator <- as.factor(mhp$elevator)
mhp$garage <- as.factor(mhp$garage)
mhp$neighborhood <- as.factor(mhp$neighborhood)
mhp$district <- as.factor(mhp$district)2.2 Tratamiento de datos faltantes
Convertimos en NA aquellas observaciones que claramente contienen algún error:
La casa de 41 habitaciones (probablemente fueran 4 y esté mal imputado)
Las que tienen menos de 3 m2 (probablemente mal imputados al usar el punto para separar los miles)
La casa con un precio de 725 (claramente equivocado).
mhp$rooms[mhp$rooms > 40] <- NA
mhp$m2[mhp$m2 < 3] <- NA
mhp$price[mhp$price < 1000] <- NAPor otro lado, la variable house_type_2 tiene 469 filas sin datos. Junto a los NA añadidos en el paso anterior son el 3.2% del total, por lo que optaremos por eliminar estas filas del dataset.
sum(is.na(mhp))/nrow(mhp)## [1] 0.03205008mhp <- na.omit(mhp)2.3 División del dataset
Para finalizar la preparación de los datos, dividiremos nuestro dataset en tres partes, una de entrenamiento (train), con el que trabajaremos durante todo el análisis, otra de prueba (test) que nos servirá para comprobar la eficacia de los diferentes modelos predictivos que se apliquen, y finalmente uno de validación (validation), que se mantendrá sin usar hasta el último día del proyecto y con el cual validaremos el modelo final.
set.seed(108)
numero_total = nrow(mhp)
# Porcentajes de train, test y validation
w_train = .5
w_test = .25
w_validation = 1 - (w_train + w_test)
# Todos los índices
indices = seq(1:numero_total)
# Muestreo
indices_train = sample(1:numero_total, numero_total * w_train)
indices_test = sample(indices[-indices_train], numero_total * w_test)
indices_validation = indices[-c(indices_train,indices_test)]
# Agrupamos
mhp_train = mhp[indices_train,]
mhp_test = mhp[indices_test,]
mhp_validation = mhp[indices_validation,]3 Análisis univariante
3.1 Análisis de variables cualitativas
Una vez preparado el dataset comenzamos con el análisis de las variables cualitativas, que son las siguientes:
Tipo de vivienda
merge(setNames(as.data.frame(table(mhp_train$house_type)), c("house_type", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$house_type))*100, 2)), c("house_type", "prop (%)"))
) %>%
arrange(desc(count)) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| house_type | count | prop (%) |
|---|---|---|
| planta 1 | 1581 | 20.45 |
| planta 2 | 1332 | 17.23 |
| planta 3 | 1180 | 15.26 |
| bajo | 996 | 12.88 |
| planta 4 | 801 | 10.36 |
| planta 5 | 545 | 7.05 |
| planta 6 | 317 | 4.10 |
| planta 7 | 196 | 2.53 |
| chalet | 180 | 2.33 |
| entreplanta | 131 | 1.69 |
| planta 8 | 119 | 1.54 |
| planta 9 | 78 | 1.01 |
| casa | 65 | 0.84 |
| semi-sotano | 48 | 0.62 |
| planta 10 | 44 | 0.57 |
| planta 12 | 23 | 0.30 |
| planta 13 | 22 | 0.28 |
| planta 11 | 21 | 0.27 |
| sotano | 19 | 0.25 |
| planta 14 | 10 | 0.13 |
| planta -1 | 8 | 0.10 |
| planta 15 | 6 | 0.08 |
| planta 16 | 4 | 0.05 |
| planta 18 | 2 | 0.03 |
| planta 19 | 2 | 0.03 |
| planta 20 | 2 | 0.03 |
ggplot(mhp_train, aes(house_type)) +
geom_bar(fill = "#0BB363") +
coord_flip() +
labs(x = "Tipo de vivienda", y = "Número de viviendas", title = "Viviendas por tipo") +
theme(plot.title = element_text(hjust = 0.5))
Como era de esperar, las viviendas más habituales son pisos en plantas no demasiado altas (bajos, primeros, segundos…). El número de casas y chalets también es pequeño en proporción.
Orientación exterior
merge(setNames(as.data.frame(table(mhp_train$exterior)), c("exterior", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$exterior))*100, 2)), c("exterior", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover"))| exterior | count | prop (%) |
|---|---|---|
| 0 | 848 | 10.97 |
| 1 | 6884 | 89.03 |
ggplot(mhp_train, aes(exterior)) +
geom_bar(fill = "#0BB363") +
labs(x = "Exterior", y = "Número de vivienda de viviendas", title = "Viviendas exteriores") +
theme(plot.title = element_text(hjust = 0.5))
Son una inmensa mayoría los pisos que dan a exterior.
Ascensor
merge(setNames(as.data.frame(table(mhp_train$elevator)), c("elevator", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$elevator))*100, 2)), c("elevator", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover"))| elevator | count | prop (%) |
|---|---|---|
| 0 | 2169 | 28.05 |
| 1 | 5563 | 71.95 |
ggplot(mhp_train, aes(elevator)) +
geom_bar(fill = "#0BB363") +
labs(x = "Ascensor", y = "Número de viviendas", title = "Viviendas con ascensor") +
theme(plot.title = element_text(hjust = 0.5))
También es más frecuente que tengan ascensor, aunque casi el 30% no lo tienen.
Garaje
merge(setNames(as.data.frame(table(mhp_train$garage)), c("garage", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$garage))*100, 2)), c("garage", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover"))| garage | count | prop (%) |
|---|---|---|
| 0 | 5570 | 72.04 |
| 1 | 2162 | 27.96 |
ggplot(mhp_train, aes(garage)) +
geom_bar(fill = "#0BB363") +
labs(x = "Garaje", y = "Número de viviendas", title = "Viviendas con garaje") +
theme(plot.title = element_text(hjust = 0.5))
Parece normal que la mayoría de pisos no tuvieran garaje incluido.
Barrio
Esta variable no será muy útil para el análisis, ya que en muchos casos son descripciones de la vivienda hechas por el usuario, por lo que habría que hacer un tratamiento previo para poder agruparlas. Algo innecesario ya que ese mismo efecto se consigue con el siguiente campo: distritos.
merge(setNames(as.data.frame(table(mhp_train$neighborhood)), c("neighborhood", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$neighborhood))*100, 2)), c("neighborhood", "prop (%)"))
) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| neighborhood | count | prop (%) |
|---|---|---|
| 12 de Octubre-Orcasur | 19 | 0.25 |
| Abrantes | 54 | 0.70 |
| Acacias | 50 | 0.65 |
| Adelfas | 30 | 0.39 |
| Aeropuerto | 2 | 0.03 |
| Águilas | 33 | 0.43 |
| Alameda de Osuna | 32 | 0.41 |
| Almagro | 140 | 1.81 |
| Almendrales | 42 | 0.54 |
| Aluche | 56 | 0.72 |
| Ambroz | 30 | 0.39 |
| Amposta | 18 | 0.23 |
| Apóstol Santiago | 19 | 0.25 |
| Arapiles | 67 | 0.87 |
| Aravaca | 47 | 0.61 |
| Arcos | 26 | 0.34 |
| Argüelles | 142 | 1.84 |
| Arroyo del Fresno | 11 | 0.14 |
| Atalaya | 6 | 0.08 |
| Ático en Almagro | 8 | 0.10 |
| Ático en Arapiles | 1 | 0.01 |
| Ático en Aravaca | 0 | 0.00 |
| Ático en Arcos | 1 | 0.01 |
| Ático en Argüelles | 0 | 0.00 |
| Ático en Bernabéu-Hispanoamérica | 1 | 0.01 |
| Ático en Campo de las Naciones-Corralejos | 1 | 0.01 |
| Ático en Castellana | 10 | 0.13 |
| Ático en Castilla | 4 | 0.05 |
| Ático en Chueca-Justicia | 1 | 0.01 |
| Ático en Ciudad Universitaria | 1 | 0.01 |
| Ático en Colina | 1 | 0.01 |
| Ático en Concepción | 1 | 0.01 |
| Ático en Conde Orgaz-Piovera | 1 | 0.01 |
| Ático en Costillares | 1 | 0.01 |
| Ático en Cuatro Caminos | 0 | 0.00 |
| Ático en Delicias | 1 | 0.01 |
| Ático en El Cañaveral | 1 | 0.01 |
| Ático en El Viso | 8 | 0.10 |
| Ático en Fuente del Berro | 0 | 0.00 |
| Ático en Gaztambide | 2 | 0.03 |
| Ático en Goya | 3 | 0.04 |
| Ático en Guindalera | 2 | 0.03 |
| Ático en Horcajo | 1 | 0.01 |
| Ático en Ibiza | 2 | 0.03 |
| Ático en Jerónimos | 1 | 0.01 |
| Ático en Las Tablas | 1 | 0.01 |
| Ático en Lista | 1 | 0.01 |
| Ático en Malasaña-Universidad | 1 | 0.01 |
| Ático en Niño Jesús | 1 | 0.01 |
| Ático en Nueva España | 4 | 0.05 |
| Ático en Nuevos Ministerios-Ríos Rosas | 3 | 0.04 |
| Ático en Pacífico | 1 | 0.01 |
| Ático en Palacio | 0 | 0.00 |
| Ático en Palomas | 1 | 0.01 |
| Ático en Peñagrande | 1 | 0.01 |
| Ático en Pinar del Rey | 0 | 0.00 |
| Ático en Pueblo Nuevo | 1 | 0.01 |
| Ático en Puerta Bonita | 0 | 0.00 |
| Ático en Puerta del Ángel | 1 | 0.01 |
| Ático en Recoletos | 3 | 0.04 |
| Ático en Rejas | 0 | 0.00 |
| Ático en Rosas | 1 | 0.01 |
| Ático en Salvador | 1 | 0.01 |
| Ático en San Juan Bautista | 1 | 0.01 |
| Ático en San Pascual | 1 | 0.01 |
| Ático en Simancas | 0 | 0.00 |
| Ático en Sol | 3 | 0.04 |
| Ático en Trafalgar | 1 | 0.01 |
| Ático en Valdeacederas | 1 | 0.01 |
| Ático en Valdebebas - Valdefuentes | 1 | 0.01 |
| Ático en Valdebernardo - Valderrivas | 0 | 0.00 |
| Ático en Valdemarín | 2 | 0.03 |
| Ático en Vallehermoso | 2 | 0.03 |
| Bellas Vistas | 89 | 1.15 |
| Bernabéu-Hispanoamérica | 91 | 1.18 |
| Berruguete | 82 | 1.06 |
| Buena Vista | 40 | 0.52 |
| Butarque | 14 | 0.18 |
| Campamento | 16 | 0.21 |
| Campo de las Naciones-Corralejos | 16 | 0.21 |
| Canillas | 43 | 0.56 |
| Canillejas | 78 | 1.01 |
| Casa de Campo | 18 | 0.23 |
| Casa o chalet independiente en Aravaca | 7 | 0.09 |
| Casa o chalet independiente en Bernabéu-Hispanoamérica | 2 | 0.03 |
| Casa o chalet independiente en Canillas | 2 | 0.03 |
| Casa o chalet independiente en Canillejas | 1 | 0.01 |
| Casa o chalet independiente en Casa de Campo | 0 | 0.00 |
| Casa o chalet independiente en Casco Histórico de Barajas | 0 | 0.00 |
| Casa o chalet independiente en Castilla | 0 | 0.00 |
| Casa o chalet independiente en Ciudad Universitaria | 5 | 0.06 |
| Casa o chalet independiente en Conde Orgaz-Piovera | 7 | 0.09 |
| Casa o chalet independiente en El Plantío | 1 | 0.01 |
| Casa o chalet independiente en El Viso | 4 | 0.05 |
| Casa o chalet independiente en Fuentelarreina | 0 | 0.00 |
| Casa o chalet independiente en Guindalera | 0 | 0.00 |
| Casa o chalet independiente en Mirasierra | 4 | 0.05 |
| Casa o chalet independiente en Niño Jesús | 1 | 0.01 |
| Casa o chalet independiente en Nueva España | 6 | 0.08 |
| Casa o chalet independiente en Opañel | 0 | 0.00 |
| Casa o chalet independiente en Orcasitas | 0 | 0.00 |
| Casa o chalet independiente en Peñagrande | 4 | 0.05 |
| Casa o chalet independiente en Rejas | 1 | 0.01 |
| Casa o chalet independiente en Tres Olivos - Valverde | 2 | 0.03 |
| Casa o chalet independiente en Valdemarín | 7 | 0.09 |
| Casa o chalet independiente en Valdezarza | 0 | 0.00 |
| Casa o chalet independiente en Villaverde Alto | 0 | 0.00 |
| Casa o chalet independiente en Vista Alegre | 0 | 0.00 |
| Casco Histórico de Barajas | 22 | 0.28 |
| Casco Histórico de Vallecas | 87 | 1.13 |
| Casco Histórico de Vicálvaro | 35 | 0.45 |
| Castellana | 172 | 2.22 |
| Castilla | 62 | 0.80 |
| Chalet adosado en Alameda de Osuna | 1 | 0.01 |
| Chalet adosado en Aravaca | 0 | 0.00 |
| Chalet adosado en Arroyo del Fresno | 1 | 0.01 |
| Chalet adosado en Bernabéu-Hispanoamérica | 1 | 0.01 |
| Chalet adosado en Buena Vista | 1 | 0.01 |
| Chalet adosado en Campo de las Naciones-Corralejos | 1 | 0.01 |
| Chalet adosado en Canillas | 1 | 0.01 |
| Chalet adosado en Canillejas | 0 | 0.00 |
| Chalet adosado en Ciudad Jardín | 1 | 0.01 |
| Chalet adosado en Ciudad Universitaria | 0 | 0.00 |
| Chalet adosado en Concepción | 1 | 0.01 |
| Chalet adosado en Conde Orgaz-Piovera | 5 | 0.06 |
| Chalet adosado en El Plantío | 2 | 0.03 |
| Chalet adosado en El Viso | 4 | 0.05 |
| Chalet adosado en Fuente del Berro | 1 | 0.01 |
| Chalet adosado en Guindalera | 0 | 0.00 |
| Chalet adosado en Los Cármenes | 0 | 0.00 |
| Chalet adosado en Mirasierra | 1 | 0.01 |
| Chalet adosado en Moscardó | 1 | 0.01 |
| Chalet adosado en Nueva España | 3 | 0.04 |
| Chalet adosado en Numancia | 0 | 0.00 |
| Chalet adosado en Orcasitas | 2 | 0.03 |
| Chalet adosado en Palomas | 1 | 0.01 |
| Chalet adosado en Peñagrande | 1 | 0.01 |
| Chalet adosado en Pinar del Rey | 1 | 0.01 |
| Chalet adosado en Puerta Bonita | 0 | 0.00 |
| Chalet adosado en Rejas | 1 | 0.01 |
| Chalet adosado en Salvador | 0 | 0.00 |
| Chalet adosado en Tres Olivos - Valverde | 0 | 0.00 |
| Chalet adosado en Valdeacederas | 0 | 0.00 |
| Chalet adosado en Valdebebas - Valdefuentes | 1 | 0.01 |
| Chalet adosado en Valdezarza | 1 | 0.01 |
| Chalet en Águilas | 0 | 0.00 |
| Chalet en Aravaca | 1 | 0.01 |
| Chalet en Ciudad Universitaria | 1 | 0.01 |
| Chalet en El Plantío | 1 | 0.01 |
| Chalet en El Viso | 0 | 0.00 |
| Chalet en Entrevías | 1 | 0.01 |
| Chalet en Fuentelarreina | 0 | 0.00 |
| Chalet en Mirasierra | 1 | 0.01 |
| Chalet en Rejas | 0 | 0.00 |
| Chalet en Valdemarín | 0 | 0.00 |
| Chalet pareado en Abrantes | 1 | 0.01 |
| Chalet pareado en Alameda de Osuna | 0 | 0.00 |
| Chalet pareado en Apóstol Santiago | 1 | 0.01 |
| Chalet pareado en Aravaca | 3 | 0.04 |
| Chalet pareado en Arroyo del Fresno | 1 | 0.01 |
| Chalet pareado en Bernabéu-Hispanoamérica | 0 | 0.00 |
| Chalet pareado en Campo de las Naciones-Corralejos | 0 | 0.00 |
| Chalet pareado en Canillas | 0 | 0.00 |
| Chalet pareado en Canillejas | 0 | 0.00 |
| Chalet pareado en Ciudad Jardín | 0 | 0.00 |
| Chalet pareado en Ciudad Universitaria | 1 | 0.01 |
| Chalet pareado en Colina | 1 | 0.01 |
| Chalet pareado en Conde Orgaz-Piovera | 4 | 0.05 |
| Chalet pareado en El Plantío | 2 | 0.03 |
| Chalet pareado en El Viso | 3 | 0.04 |
| Chalet pareado en Fuentelarreina | 0 | 0.00 |
| Chalet pareado en Guindalera | 0 | 0.00 |
| Chalet pareado en Mirasierra | 3 | 0.04 |
| Chalet pareado en Nueva España | 0 | 0.00 |
| Chalet pareado en Palomas | 4 | 0.05 |
| Chalet pareado en Peñagrande | 1 | 0.01 |
| Chalet pareado en Rejas | 1 | 0.01 |
| Chalet pareado en Tres Olivos - Valverde | 0 | 0.00 |
| Chalet pareado en Valdebebas - Valdefuentes | 1 | 0.01 |
| Chopera | 34 | 0.44 |
| Chueca-Justicia | 67 | 0.87 |
| Ciudad Jardín | 47 | 0.61 |
| Ciudad Universitaria | 40 | 0.52 |
| Colina | 25 | 0.32 |
| Comillas | 52 | 0.67 |
| Concepción | 49 | 0.63 |
| Conde Orgaz-Piovera | 35 | 0.45 |
| Costillares | 53 | 0.69 |
| Cuatro Caminos | 98 | 1.27 |
| Cuatro Vientos | 1 | 0.01 |
| Cuzco-Castillejos | 98 | 1.27 |
| Delicias | 58 | 0.75 |
| Dúplex en Almagro | 6 | 0.08 |
| Dúplex en Apóstol Santiago | 0 | 0.00 |
| Dúplex en Arapiles | 0 | 0.00 |
| Dúplex en Aravaca | 3 | 0.04 |
| Dúplex en Arcos | 0 | 0.00 |
| Dúplex en Argüelles | 0 | 0.00 |
| Dúplex en Arroyo del Fresno | 1 | 0.01 |
| Dúplex en Bernabéu-Hispanoamérica | 0 | 0.00 |
| Dúplex en Berruguete | 1 | 0.01 |
| Dúplex en Butarque | 0 | 0.00 |
| Dúplex en Canillas | 1 | 0.01 |
| Dúplex en Canillejas | 1 | 0.01 |
| Dúplex en Casco Histórico de Barajas | 0 | 0.00 |
| Dúplex en Casco Histórico de Vallecas | 0 | 0.00 |
| Dúplex en Castellana | 1 | 0.01 |
| Dúplex en Castilla | 1 | 0.01 |
| Dúplex en Chueca-Justicia | 1 | 0.01 |
| Dúplex en Colina | 1 | 0.01 |
| Dúplex en Concepción | 2 | 0.03 |
| Dúplex en Conde Orgaz-Piovera | 1 | 0.01 |
| Dúplex en Costillares | 0 | 0.00 |
| Dúplex en Cuzco-Castillejos | 1 | 0.01 |
| Dúplex en El Viso | 4 | 0.05 |
| Dúplex en Entrevías | 1 | 0.01 |
| Dúplex en Fontarrón | 1 | 0.01 |
| Dúplex en Fuente del Berro | 0 | 0.00 |
| Dúplex en Goya | 1 | 0.01 |
| Dúplex en Imperial | 1 | 0.01 |
| Dúplex en Jerónimos | 4 | 0.05 |
| Dúplex en Las Tablas | 2 | 0.03 |
| Dúplex en Legazpi | 2 | 0.03 |
| Dúplex en Lista | 0 | 0.00 |
| Dúplex en Malasaña-Universidad | 0 | 0.00 |
| Dúplex en Mirasierra | 1 | 0.01 |
| Dúplex en Moscardó | 1 | 0.01 |
| Dúplex en Nueva España | 8 | 0.10 |
| Dúplex en Nuevos Ministerios-Ríos Rosas | 1 | 0.01 |
| Dúplex en Numancia | 2 | 0.03 |
| Dúplex en Opañel | 1 | 0.01 |
| Dúplex en Pacífico | 3 | 0.04 |
| Dúplex en Palacio | 1 | 0.01 |
| Dúplex en Pinar del Rey | 0 | 0.00 |
| Dúplex en Prosperidad | 1 | 0.01 |
| Dúplex en Puerta Bonita | 0 | 0.00 |
| Dúplex en Puerta del Ángel | 0 | 0.00 |
| Dúplex en Quintana | 0 | 0.00 |
| Dúplex en Recoletos | 3 | 0.04 |
| Dúplex en Rejas | 0 | 0.00 |
| Dúplex en San Isidro | 1 | 0.01 |
| Dúplex en San Juan Bautista | 6 | 0.08 |
| Dúplex en Sanchinarro | 1 | 0.01 |
| Dúplex en Simancas | 2 | 0.03 |
| Dúplex en Trafalgar | 2 | 0.03 |
| Dúplex en Tres Olivos - Valverde | 1 | 0.01 |
| Dúplex en Valdeacederas | 3 | 0.04 |
| Dúplex en Valdebernardo - Valderrivas | 0 | 0.00 |
| Dúplex en Valdemarín | 4 | 0.05 |
| Dúplex en Vallehermoso | 1 | 0.01 |
| Dúplex en Ventilla-Almenara | 0 | 0.00 |
| Dúplex en Villaverde Alto | 1 | 0.01 |
| El Cañaveral | 50 | 0.65 |
| El Pardo | 3 | 0.04 |
| El Plantío | 3 | 0.04 |
| El Viso | 109 | 1.41 |
| Ensanche de Vallecas - La Gavia | 60 | 0.78 |
| Entrevías | 38 | 0.49 |
| Estrella | 27 | 0.35 |
| Fontarrón | 31 | 0.40 |
| Fuente del Berro | 45 | 0.58 |
| Fuentelarreina | 9 | 0.12 |
| Gaztambide | 77 | 1.00 |
| Goya | 257 | 3.32 |
| Guindalera | 97 | 1.25 |
| Hellín | 17 | 0.22 |
| Horcajo | 3 | 0.04 |
| Huertas-Cortes | 48 | 0.62 |
| Ibiza | 56 | 0.72 |
| Imperial | 48 | 0.62 |
| Jerónimos | 62 | 0.80 |
| La Paz | 48 | 0.62 |
| Las Tablas | 39 | 0.50 |
| Lavapiés-Embajadores | 115 | 1.49 |
| Legazpi | 38 | 0.49 |
| Lista | 134 | 1.73 |
| Los Ángeles | 43 | 0.56 |
| Los Berrocales | 3 | 0.04 |
| Los Cármenes | 50 | 0.65 |
| Los Rosales | 34 | 0.44 |
| Lucero | 43 | 0.56 |
| Malasaña-Universidad | 126 | 1.63 |
| Marroquina | 17 | 0.22 |
| Media Legua | 20 | 0.26 |
| Mirasierra | 15 | 0.19 |
| Montecarmelo | 14 | 0.18 |
| Moscardó | 47 | 0.61 |
| Niño Jesús | 24 | 0.31 |
| Nueva España | 76 | 0.98 |
| Nuevos Ministerios-Ríos Rosas | 91 | 1.18 |
| Numancia | 90 | 1.16 |
| Opañel | 54 | 0.70 |
| Orcasitas | 38 | 0.49 |
| Pacífico | 83 | 1.07 |
| Palacio | 88 | 1.14 |
| Palomas | 20 | 0.26 |
| Palomeras Bajas | 42 | 0.54 |
| Palomeras sureste | 37 | 0.48 |
| Palos de Moguer | 73 | 0.94 |
| Pau de Carabanchel | 26 | 0.34 |
| Pavones | 6 | 0.08 |
| Peñagrande | 77 | 1.00 |
| Pilar | 63 | 0.81 |
| Pinar del Rey | 67 | 0.87 |
| Portazgo | 43 | 0.56 |
| Pradolongo | 32 | 0.41 |
| Prosperidad | 86 | 1.11 |
| Pueblo Nuevo | 132 | 1.71 |
| Puerta Bonita | 89 | 1.15 |
| Puerta del Ángel | 137 | 1.77 |
| Quintana | 68 | 0.88 |
| Recoletos | 181 | 2.34 |
| Rejas | 59 | 0.76 |
| Rosas | 30 | 0.39 |
| Salvador | 21 | 0.27 |
| San Cristóbal | 30 | 0.39 |
| San Diego | 120 | 1.55 |
| San Fermín | 33 | 0.43 |
| San Isidro | 91 | 1.18 |
| San Juan Bautista | 39 | 0.50 |
| San Pascual | 48 | 0.62 |
| Sanchinarro | 61 | 0.79 |
| Santa Eugenia | 16 | 0.21 |
| Simancas | 60 | 0.78 |
| Sol | 66 | 0.85 |
| Timón | 14 | 0.18 |
| Trafalgar | 97 | 1.25 |
| Tres Olivos - Valverde | 51 | 0.66 |
| Valdeacederas | 104 | 1.35 |
| Valdebebas - Valdefuentes | 42 | 0.54 |
| Valdebernardo - Valderrivas | 22 | 0.28 |
| Valdemarín | 19 | 0.25 |
| Valdezarza | 42 | 0.54 |
| Vallehermoso | 72 | 0.93 |
| Ventas | 103 | 1.33 |
| Ventilla-Almenara | 46 | 0.59 |
| Villaverde Alto | 87 | 1.13 |
| Vinateros | 20 | 0.26 |
| Virgen del Cortijo - Manoteras | 31 | 0.40 |
| Vista Alegre | 102 | 1.32 |
| Zofío | 25 | 0.32 |
Distrito
En este caso los datos si están correctamente recogidos, repartiendo las casas entre los 21 distritos de Madrid.Casi todos con una representación considerable, en lo que destaca el barrio de Salamanca, con casi el doble de viviendas que el segundo que más tiene, Chamberí.
merge(setNames(as.data.frame(table(mhp_train$district)), c("district", "count")),
setNames(as.data.frame(round(prop.table(table(mhp_train$district))*100, 2)), c("district", "prop (%)"))
) %>%
arrange(desc(count)) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| district | count | prop (%) |
|---|---|---|
| barrio de salamanca | 911 | 11.78 |
| chamberi | 571 | 7.38 |
| ciudad lineal | 540 | 6.98 |
| chamartin | 526 | 6.80 |
| tetuan | 523 | 6.76 |
| centro | 517 | 6.69 |
| carabanchel | 512 | 6.62 |
| puente-de-vallecas | 374 | 4.84 |
| fuencarral | 356 | 4.60 |
| moncloa | 353 | 4.57 |
| hortaleza | 352 | 4.55 |
| latina | 337 | 4.36 |
| san-blas | 319 | 4.13 |
| Arganzuela | 305 | 3.94 |
| retiro | 295 | 3.82 |
| usera | 240 | 3.10 |
| villaverde | 209 | 2.70 |
| villa-de-vallecas | 163 | 2.11 |
| vicalvaro | 141 | 1.82 |
| moratalaz | 99 | 1.28 |
| barajas | 89 | 1.15 |
ggplot(mhp_train, aes(district)) +
geom_bar(fill = "#0BB363") +
labs(x = "Distrito", y = "Número de viviendas", title = "Viviendas por distrito") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
3.2 Análisis de variables cuantitativas
Precio
data.frame(summarise(mhp_train,
min = min(price),
max = max(price),
median = median(price),
mean = mean(price),
sd = sd(price))) %>%
kbl() %>%
kable_material(c("striped", "hover"))| min | max | median | mean | sd |
|---|---|---|---|---|
| 47500 | 8950000 | 369000 | 631920.1 | 753136.6 |
ggplot(mhp_train, aes(y = price)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Precio", title = "Boxplot de precios") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(price)) +
geom_histogram(aes(y = ..count..), bins = 50, position = "dodge", fill = "#0BB363") +
labs(x = "Precio", y = "Número de viviendas", title = "Viviendas por Precio") +
scale_x_continuous(breaks = pretty(mhp_train$price, n = 10), labels = comma_format()) +
scale_y_continuous(labels = function(x) sprintf("%.0f", x)) +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30))
Habitaciones
data.frame(summarise(mhp_train,
min = min(rooms),
max = max(rooms),
median = median(rooms),
mean = mean(rooms),
sd = sd(rooms))) %>%
kbl() %>%
kable_material(c("striped", "hover"))| min | max | median | mean | sd |
|---|---|---|---|---|
| 1 | 17 | 3 | 2.848552 | 1.302995 |
ggplot(mhp_train, aes(y = rooms)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Habitaciones", title = "Boxplot de habitaciones") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(rooms)) +
geom_histogram(aes(), bins = 16, position = "dodge", fill = "#0BB363") +
geom_density(alpha=.2, fill = "red") +
labs(x = "Habitaciones", y = "Número de viviendas", title = "Viviendas por número de habitaciones") +
theme(plot.title = element_text(hjust = 0.5))
m²
data.frame(summarise(mhp_train,
min = min(m2),
max = max(m2),
median = median(m2),
mean = mean(m2),
sd = sd(m2))) %>%
kbl() %>%
kable_material(c("striped", "hover"))| min | max | median | mean | sd |
|---|---|---|---|---|
| 20 | 986 | 95 | 127.1208 | 102.497 |
ggplot(mhp_train, aes(y = m2)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Metros cuadrados", title = "Boxplot de superficies") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(m2)) +
geom_histogram(aes(y=..count..), bins = 50, position = "dodge", fill = "#0BB363") +
geom_density(alpha=.3, fill = "red") +
labs(x = "Metros cuadrados", y = "Número de viviendas", title = "Viviendas por metros cuadrados") +
theme(plot.title = element_text(hjust = 0.5))
4 Análisis multivariante
Una vez analizadas todas las variables de forma individual, podemos buscar algunas relaciones entre ellas.
Como se ve en los gráficos inferiores hay una relción fuerte entre el precio, los metros cuadradados y el número de habitaciones.
plot_grid(
ggcorrplot(cor(mhp_train %>% select_if(is.numeric)),type = "lower", lab=TRUE),
ggplot(mhp_train, aes(x = m2, y = price)) +
geom_point() +
geom_smooth() +
ggtitle('Reación precio y m²') +
theme(plot.title = element_text(hjust = 0.5)),
ggplot(mhp_train, aes(x = rooms, y = price)) +
geom_point() +
geom_smooth() +
ggtitle('Reación precio y habitaciones') +
theme(plot.title = element_text(hjust = 0.5)),
ggplot(mhp_train, aes(x = rooms, y = m2)) +
geom_point() +
geom_smooth() +
ggtitle('Reación m² y habitaciones') +
theme(plot.title = element_text(hjust = 0.5)),
nrow = 2
)
Introduciendo algunas de las variables categóricas se puede observar como varía el precio con la superficie o el número de habitaciones dependiendo de si incluye garaje, ascensor u orientación exterior.
plot_grid(
ggplot(mhp_train, aes(x = m2, y = price, colour = exterior)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-superficie, exterior') ,
ggplot(mhp_train, aes(x = rooms, y = price, colour = exterior)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-habitaciones, exterior') ,
nrow =1
)
plot_grid(
ggplot(mhp_train, aes(x = m2, y = price, colour = elevator)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-superficie, ascensor'),
ggplot(mhp_train, aes(x = rooms, y = price, colour = elevator)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-habitación, ascensor'),
nrow =1
)
plot_grid(
ggplot(mhp_train, aes(x = m2, y = price, colour = garage)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-superficie, garaje'),
ggplot(mhp_train, aes(x = rooms, y = price, colour = garage)) +
geom_point() +
geom_smooth(method = "lm") +
scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
ggtitle('Rel. precio-habitaciones, garaje'),
nrow =1
)
Representando boxplot por distrito se ve claramente que los barrios ricos como Salamanca o Chamberí, a parte de tener precios medios más altos, tienen mayor número de viviendas con precios atípicos (por encima), mientras que en el número de habitaciones no se aprecia una diferencia tan notoria.
ggplot(mhp_train, aes(district, price)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Precio", title = "Boxplot de precios por distrito") +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(district, rooms)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Habitaciones", title = "Boxplot de habitaciones por distrito") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
ggplot(mhp_train, aes(district, m2)) +
geom_boxplot(fill = "#0BB363") +
labs(y = "Metros cuadrados", title = "Boxplot de superficie por distrito") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))
En las siguientes figuras se representa el precio medio por distrito y como varía con respecto a otras variables de interés.
mhp_train %>%
group_by(district, m2) %>%
summarize(avg_price = mean(price)) %>%
ggplot(aes(x = m2, y = avg_price)) +
geom_point(size = 0.5) +
facet_wrap(~ district) +
scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))
mhp_train %>%
group_by(district, exterior) %>%
summarise(avg_price = mean(price)) %>%
ggplot(aes(x=district, y=avg_price, fill=exterior)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
ggtitle("Precio medio por distrito y exterior") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(x = "Distrito", y = "Precio medio") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
mhp_train %>%
group_by(district, elevator) %>%
summarise(avg_price = mean(price)) %>%
ggplot(aes(x=district, y=avg_price, fill=elevator)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
ggtitle("Precio medio por distrito y ascensor") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(x = "Distrito", y = "Precio medio") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
mhp_train %>%
group_by(district, garage) %>%
summarise(avg_price = mean(price)) %>%
ggplot(aes(x=district, y=avg_price, fill=garage)) +
geom_bar(stat = "identity", position = "dodge") +
scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
ggtitle("Precio medio por distrito y garaje") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(x = "Distrito", y = "Precio medio") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
5 Procesado de variables cualitativas
mhp_train %>%
group_by(house_type) %>%
summarise(precio_medio = mean(price)) %>%
arrange(precio_medio)## # A tibble: 26 × 2
## house_type precio_medio
## <fct> <dbl>
## 1 sotano 279031.
## 2 planta -1 298188.
## 3 semi-sotano 302916.
## 4 bajo 316306.
## 5 planta 11 458019.
## 6 planta 1 540375.
## 7 planta 2 595355.
## 8 entreplanta 598834.
## 9 planta 10 605368.
## 10 planta 4 605711.
## # … with 16 more rowsPara el posterior análisis transformamos la columna house_type, agrupando en “piso” (independientemente de la planta en la que esté), “casa” (casa o chalet) y “sotano”.
mhp_train <- mhp_train %>%
mutate(tipo_casa = case_when(
grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
TRUE ~ "piso"
))
mhp_train$tipo_casa <- as.factor(mhp_train$tipo_casa)También se agruparán los distritos en función del precio medio del mismo: caro, medio o barato.
mhp_train <- mhp_train %>%
group_by(district) %>%
mutate(precio_distrito = mean(price)) %>%
ungroup()
# Calcular los percentiles
percentiles <- quantile(mhp_train$precio_distrito, probs = c(0, 1/3, 2/3, 1))
# Asignar etiquetas
mhp_train$precio_distrito <- cut(mhp_train$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)Además, de cara a desarrollar algunos modelos, se creará una variable target binaria que dividirá las casas entre “cara” y “barata” (1/0).
mhp_train$precio_bin <- ifelse(mhp_train$price > median(mhp_train$price), "cara", "barata")
mhp_train$precio_bin = ifelse(mhp_train$precio_bin == "cara", 1, 0)
mhp_train$precio_bin = factor(mhp_train$precio_bin, levels = c(1,0))# Aplicamos todos los cambios también a test.
# Categorizamos house_type
mhp_test <- mhp_test %>%
mutate(tipo_casa = case_when(
grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
TRUE ~ "piso"
))
mhp_test$tipo_casa <- as.factor(mhp_test$tipo_casa)
# Categorizamos los distritos
mhp_test <- mhp_test %>%
group_by(district) %>%
mutate(precio_distrito = mean(price)) %>%
ungroup()
percentiles <- quantile(mhp_test$precio_distrito, probs = c(0, 1/3, 2/3, 1))
mhp_test$precio_distrito <- cut(mhp_test$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)
# Creamos variable target
mhp_test$precio_bin <- ifelse(mhp_test$price > median(mhp_test$price), "cara", "barata")
mhp_test$precio_bin = ifelse(mhp_test$precio_bin == "cara", 1, 0)
mhp_test$precio_bin = factor(mhp_test$precio_bin, levels = c(1,0))6 Modelo de regresión lineal
Con el dataset ya preparado, se crea un modelo de regresión lineal.
En este caso, tras hacer diversas pruebas se ha optado por introducir sólo 3 variables (tipo_casa, m2 y precio_distrito), alcanzando un R2 de 0.7553.
Esta cifra podía incrementarse algo, sin embargo es a costa de introducir más variables y por tanto más complejidad al modelo, haciéndolo menos explicable.
lm_fit <- lm(price ~ m2+tipo_casa+precio_distrito, data=mhp_train)
summary(lm_fit)##
## Call:
## lm(formula = price ~ m2 + tipo_casa + precio_distrito, data = mhp_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3277290 -137875 -17495 95705 5483233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -915949.86 30574.68 -29.958 <2e-16 ***
## m2 6052.79 51.28 118.027 <2e-16 ***
## tipo_casapiso 673879.37 27865.36 24.183 <2e-16 ***
## tipo_casasotano 606283.79 28773.95 21.071 <2e-16 ***
## precio_distritomedio 101183.27 10337.15 9.788 <2e-16 ***
## precio_distritocaro 373747.30 11419.31 32.729 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 372600 on 7726 degrees of freedom
## Multiple R-squared: 0.7553, Adjusted R-squared: 0.7552
## F-statistic: 4771 on 5 and 7726 DF, p-value: < 2.2e-16coef(lm_fit)## (Intercept) m2 tipo_casapiso
## -915949.864 6052.795 673879.371
## tipo_casasotano precio_distritomedio precio_distritocaro
## 606283.787 101183.273 373747.298residuals=lm_fit$residuals
autoplot(lm_fit)
Observando los residuos se ve como hay cierta heterocedasticidad y en la gráfica QQ no se ajustan a la normal, por lo que el módelo no es bueno.
# Realizar predicciones en la partición test
predictions <- predict(lm_fit, newdata = mhp_test)
# Comparar las predicciones con los valores reales
rmse <- rmse(predictions, mhp_test$price)
mae <- mae(predictions, mhp_test$price)
r2 <- R2(predictions, mhp_test$price)
# Imprimir las métricas de evaluación
cat(paste0("RMSE: ", round(rmse, 2), "\n"))## RMSE: 398263.79cat(paste0("MAE: ", round(mae, 2), "\n"))## MAE: 199633.64cat(paste0("R-squared: ", round(r2, 2)))## R-squared: 0.73Un R2 de 0.73 (bastante similar al calculado en train) indica que se explica con esa probabilidad la variabilidad de la variable dependiente, por lo tanto podríamos considerarlo razonablemente útil, sin embargo sería preferible mejorar los valores del error de predicción promedio o el cuadrático.
7 Conclusiones preliminares
El análisis realizado hasta el momento ha servido para familiarizarse con el mercado inmobiliario de Madrid (o al menos parte de él) y conocer como se comportan algunas de las variables más relevantes en el valor de las viviendas.
Con estos conocimientos se abordarán nuevos modelos para intentar conseguir mejores predicciones y simultaneamente comprender y desarrollar estos mismos modelos.
str(mhp_train)## tibble [7,732 × 12] (S3: tbl_df/tbl/data.frame)
## $ price : num [1:7732] 499000 385000 950000 460000 335000 580000 940000 85300 399000 125000 ...
## $ house_type : Factor w/ 26 levels "bajo","casa",..: 21 16 20 10 18 19 16 16 6 1 ...
## $ exterior : Factor w/ 2 levels "1","0": 1 1 1 1 1 1 1 1 1 1 ...
## $ rooms : num [1:7732] 4 1 3 3 4 3 2 3 3 2 ...
## $ m2 : num [1:7732] 131 60 160 114 122 120 124 75 86 62 ...
## $ elevator : Factor w/ 2 levels "1","0": 1 2 1 1 1 1 1 1 2 2 ...
## $ garage : Factor w/ 2 levels "1","0": 2 2 2 1 2 2 2 2 2 2 ...
## $ neighborhood : Factor w/ 341 levels "12 de Octubre-Orcasur",..: 76 17 8 284 318 265 264 316 192 11 ...
## $ district : Factor w/ 21 levels "Arganzuela","barajas",..: 17 12 7 13 18 3 3 21 1 19 ...
## $ tipo_casa : Factor w/ 3 levels "casa","piso",..: 2 2 2 2 2 2 2 2 3 3 ...
## $ precio_distrito: Factor w/ 3 levels "barato","medio",..: 2 3 3 1 1 3 3 1 1 1 ...
## $ precio_bin : Factor w/ 2 levels "1","0": 1 1 1 1 2 1 1 2 1 2 ...8 Aprendizaje no supervisado
8.1 K-Means
El algoritmo K-Means es un método de aprendizaje automático no supervisado utilizado para el análisis de agrupamiento o “clustering”. El objetivo principal de este algoritmo es dividir los datos en K grupos (clusters), de modo que los puntos dentro de un mismo grupo sean lo más similares posible.
Por como funciona este método, es habitual utilizarlo en la etapa de EDA para agrupar y etiquetar observaciones con un criterior más científico que el basado en el conocimiento personal del dominio, más que para hacer predicciones (aunque también puede hacerse).
En nuestro lo utilizaremos para comprobar si hay alguna forma mejor de agrupar las casas de como lo hemos hecho anteriormente.
mhp_train_kmeans <- subset(mhp_train, select = c("m2", "rooms", "elevator", "garage", "exterior"))
mhp_train_kmeans$exterior <- as.numeric(mhp_train_kmeans$exterior)
mhp_train_kmeans$elevator <- as.numeric(mhp_train_kmeans$elevator)
mhp_train_kmeans$garage <- as.numeric(mhp_train_kmeans$garage)
# Normalizamos los datos
mhp_train_kmeans <- scale(mhp_train_kmeans)Tras normalizar los datos, aplicamos k-means sólo para las variables numéricas y mediante los siguientes test decidiremos cual es el número óptimo de clusters.
fviz_nbclust(mhp_train_kmeans, kmeans, method = "wss")
fviz_nbclust(mhp_train_kmeans, kmeans, method = "silhouette")
Utilizando la regla del codo en el gráfico wss podríamos decir que está en el 6, algo que confirma el segundo gráfico. Por lo tanto, aplicamos el modelo con K = 6
# Definir el número de grupos
k <- 6
# Ejecutar k-means
grupos <- kmeans(mhp_train_kmeans, k)
mhp_train$cluster <- as.factor(grupos$cluster)
# Graficar la relación entre las características y los grupos
ggplot(mhp_train, aes(x=m2, y=price, color=cluster)) +
geom_point() +
ggtitle("Grupos de casas en venta") +
theme(plot.title = element_text(hjust = 0.5)) +
xlab("Tamaño en m2") +
ylab("Precio")
fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "euclid",repel = TRUE,star.plot = TRUE)
fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "norm")
Graficamente es difícil diferenciar las agrupaciones, por lo que en la tabla inferior obtenemos los precios medios de cada cluster, donde se puede apreciar que la gran mayoría de agrupaciones son coherentes (salvo tal vez el 2 y el 3 que son muy similares).
También se muestra el tipo de casa y el distrito predominante en cada cluster.
# Precio medio por cluster
mhp_train %>%
group_by(cluster) %>%
summarise(precio_medio = mean(price))## # A tibble: 6 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 184203.
## 2 2 640443.
## 3 3 617779.
## 4 4 223476.
## 5 5 2228444.
## 6 6 355840.# Tipo de casa predominante por cluster
mhp_train %>%
group_by(cluster, tipo_casa) %>%
summarise(frecuencia = n()) %>%
arrange(cluster, desc(frecuencia)) %>%
slice(1) %>%
select(cluster, tipo_casa_predominante = tipo_casa)## # A tibble: 6 × 2
## # Groups: cluster [6]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 sotano
## 2 2 piso
## 3 3 piso
## 4 4 piso
## 5 5 piso
## 6 6 piso# Distrito predominante por cluster
mhp_train %>%
group_by(cluster, district) %>%
summarise(frecuencia = n()) %>%
arrange(cluster, desc(frecuencia)) %>%
slice(1) %>%
select(cluster, distrito_predominante = district)## # A tibble: 6 × 2
## # Groups: cluster [6]
## cluster distrito_predominante
## <fct> <fct>
## 1 1 tetuan
## 2 2 hortaleza
## 3 3 barrio de salamanca
## 4 4 carabanchel
## 5 5 barrio de salamanca
## 6 6 barrio de salamancaSabemos que tanto la regla del codo como cualquier otro método para definir el número de clusters no siempre es adecuada, por lo que a continuación se prueba con un bucle que prueba con número de clusters menores que 6 y muestra las agrupaciones por precio y tipo de piso anteriores.
# Ejecuta k-means para diferentes valores de K
for (k in 2:6) {
# Ejecutar k-means
grupos <- kmeans(mhp_train_kmeans, centers = k)
mhp_train$cluster <- as.factor(grupos$cluster)
print(k)
print(mhp_train %>%
group_by(cluster) %>%
summarise(precio_medio = mean(price)))
print(mhp_train %>%
group_by(cluster, tipo_casa) %>%
summarise(frecuencia = n()) %>%
arrange(cluster, desc(frecuencia)) %>%
slice(1) %>%
select(cluster, tipo_casa_predominante = tipo_casa))
}## [1] 2
## # A tibble: 2 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 394442.
## 2 2 1094473.
## # A tibble: 2 × 2
## # Groups: cluster [2]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## [1] 3
## # A tibble: 3 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 402986.
## 2 2 605000.
## 3 3 1941976.
## # A tibble: 3 × 2
## # Groups: cluster [3]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 piso
## [1] 4
## # A tibble: 4 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 579636.
## 2 2 352563.
## 3 3 665255.
## 4 4 189057.
## # A tibble: 4 × 2
## # Groups: cluster [4]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 piso
## 4 4 sotano
## [1] 5
## # A tibble: 5 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 435798.
## 2 2 241929.
## 3 3 189057.
## 4 4 665853.
## 5 5 1005091.
## # A tibble: 5 × 2
## # Groups: cluster [5]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 sotano
## 4 4 piso
## 5 5 piso
## [1] 6
## # A tibble: 6 × 2
## cluster precio_medio
## <fct> <dbl>
## 1 1 1179261.
## 2 2 652498.
## 3 3 355840.
## 4 4 217255.
## 5 5 2356279.
## 6 6 451222.
## # A tibble: 6 × 2
## # Groups: cluster [6]
## cluster tipo_casa_predominante
## <fct> <fct>
## 1 1 piso
## 2 2 piso
## 3 3 piso
## 4 4 piso
## 5 5 piso
## 6 6 pisoObservando las diferencias entre los diferentes modelos, el número de clusters que mejor diferencia las casas por su precio es 3, algo que encaja con la presunción que inicial que hicimos para dividir los tipos de casa en 3 tipos.
Si observamos además el tipo de casas más común por cluster, casi siempre es piso (salvo en los casos de clusters con muy poco precio medio, que es sótano), lo cual aunque tiene sentido porque este tipo es claramente mayoritario en todo el dataset, no nos permite apreciar a simple vista donde se acumulan el resto de viviendas.
Por lo tanto, entre los datos extraídos del análisis de clusters y el conocimiento previo que teníamos del dominio, tiene sentido mantener los 3 tipos de viviendas principales como variable para los siguientes modelos (sótano, piso y casa).
9 Técnicas de reducción de la dimensionalidad
9.1 PCA
El Análisis de Componentes Principales es una técnica de reducción de la dimensionalidad que se utiliza para transformar un conjunto de datos de múltiples variables en un nuevo conjunto de datos con un número menor de variables cuyo principal objetivo es simplificar la representación de datos mientras se conserva la mayor cantidad de información y variabilidad posible identificando patrones y estructuras subyacentes, reduciendo la multicolinealidad y el sobreajuste.
PAra nuestro dataset, seleccionaremos las columnas numéricas, descartando la columna precio puesto que será de alguna manera nuestra variable dependiente. Además creamos 3 variables dummy usando variables categóricas de tipo 0|1.
train_pca <- mhp_train[,c(3:7,10,11,12)]dummy_garaje <- model.matrix(~0 + garage,data = train_pca)
dummy_exterior <- model.matrix(~0 + exterior, data = train_pca)
dummy_ascensor <- model.matrix(~0 + elevator, data = train_pca)
dummy_tipo_casa <- model.matrix(~0 + tipo_casa, data = train_pca)
dummy_distrito <- model.matrix(~0 + precio_distrito, data = train_pca)dummy_data <- cbind(train_pca[, c("rooms", "m2")],dummy_exterior,dummy_garaje,dummy_ascensor,dummy_distrito,dummy_tipo_casa)prcomp(dummy_data)## Standard deviations (1, .., p=14):
## [1] 1.025015e+02 9.457105e-01 7.267966e-01 6.257613e-01 6.127629e-01
## [6] 5.782676e-01 4.571490e-01 4.037562e-01 1.716715e-01 2.268386e-15
## [11] 1.853843e-15 1.003401e-15 8.037840e-16 4.438592e-16
##
## Rotation (n x k) = (14 x 14):
## PC1 PC2 PC3 PC4
## rooms -8.797027e-03 -0.9916184565 -0.014616928 0.094846395
## m2 -9.999557e-01 0.0085416356 0.001971791 -0.001450128
## exterior1 -5.771447e-04 -0.0453751425 -0.009449951 -0.063478191
## exterior0 5.771447e-04 0.0453751425 0.009449951 0.063478191
## garage1 -1.538147e-03 0.0421269404 -0.223561300 0.037355645
## garage0 1.538147e-03 -0.0421269404 0.223561300 -0.037355645
## elevator1 -4.274374e-04 0.0241411854 -0.476833313 0.342257253
## elevator0 4.274374e-04 -0.0241411854 0.476833313 -0.342257253
## precio_distritobarato 1.400700e-03 -0.0385563529 0.200000119 -0.359200946
## precio_distritomedio 6.834450e-05 -0.0051663464 -0.051807433 0.201186914
## precio_distritocaro -1.469044e-03 0.0437226993 -0.148192687 0.158014031
## tipo_casacasa -8.085078e-04 0.0005892036 0.067961955 -0.037840260
## tipo_casapiso 2.818685e-05 -0.0469938502 -0.466029098 -0.505128722
## tipo_casasotano 7.803209e-04 0.0464046466 0.398067143 0.542968982
## PC5 PC6 PC7 PC8
## rooms -0.0132115607 -0.0291253564 0.074126853 -0.029792435
## m2 0.0009596296 0.0005665667 -0.002313108 -0.001386419
## exterior1 -0.1552630611 -0.0333936354 -0.327138796 0.601052285
## exterior0 0.1552630611 0.0333936354 0.327138796 -0.601052285
## garage1 -0.5627134868 -0.1309300044 0.333444368 0.033280690
## garage0 0.5627134868 0.1309300044 -0.333444368 -0.033280690
## elevator1 0.0482119831 -0.1533456743 -0.318624110 -0.141577832
## elevator0 -0.0482119831 0.1533456743 0.318624110 0.141577832
## precio_distritobarato -0.2366848900 -0.4502445095 -0.365221920 -0.322493702
## precio_distritomedio -0.1969164022 0.7585991702 -0.088292084 -0.036423761
## precio_distritocaro 0.4336012922 -0.3083546607 0.453514004 0.358917463
## tipo_casacasa -0.0381883034 0.0229829599 0.061657352 -0.001868362
## tipo_casapiso 0.1388358808 0.1326235059 0.004437444 -0.006913037
## tipo_casasotano -0.1006475774 -0.1556064657 -0.066094796 0.008781399
## PC9 PC10 PC11 PC12
## rooms 0.000862928 5.877890e-16 1.174221e-15 7.345674e-17
## m2 -0.001036478 -4.349751e-20 3.574359e-18 -4.762673e-18
## exterior1 0.017792902 8.077407e-03 2.416838e-02 -9.348594e-03
## exterior0 -0.017792902 8.077407e-03 2.416838e-02 -9.348594e-03
## garage1 -0.027676182 6.464949e-03 5.102210e-03 -8.893342e-02
## garage0 0.027676182 6.464949e-03 5.102210e-03 -8.893342e-02
## elevator1 0.086614598 6.353245e-01 -3.028358e-01 -6.708710e-02
## elevator0 -0.086614598 6.353245e-01 -3.028358e-01 -6.708710e-02
## precio_distritobarato -0.004863213 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritomedio -0.010432322 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritocaro 0.015295535 -9.943406e-02 -8.213108e-02 -5.584679e-01
## tipo_casacasa 0.809211511 -2.329879e-01 -5.148217e-01 1.145226e-01
## tipo_casapiso -0.393939361 -2.329879e-01 -5.148217e-01 1.145226e-01
## tipo_casasotano -0.415272150 -2.329879e-01 -5.148217e-01 1.145226e-01
## PC13 PC14
## rooms -9.651041e-17 -5.224615e-18
## m2 4.139246e-18 -4.009145e-19
## exterior1 5.278108e-02 -7.046115e-01
## exterior0 5.278108e-02 -7.046115e-01
## garage1 -6.995888e-01 -5.097578e-02
## garage0 -6.995888e-01 -5.097578e-02
## elevator1 1.228500e-02 -1.293890e-03
## elevator0 1.228500e-02 -1.293890e-03
## precio_distritobarato 6.884856e-02 8.609917e-03
## precio_distritomedio 6.884856e-02 8.609917e-03
## precio_distritocaro 6.884856e-02 8.609917e-03
## tipo_casacasa -1.877163e-02 -2.325502e-02
## tipo_casapiso -1.877163e-02 -2.325502e-02
## tipo_casasotano -1.877163e-02 -2.325502e-02Podemos observar que en la PC1 lleva un peso casi de 100 la variable m2, mientras que en PC2 pasa lo mismo con la variable rooms. En el resto de PCs el análisis de componentes principales está algo más repartido.
summary(prcomp(dummy_data))## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 102.5015 0.94571 0.72680 0.62576 0.61276 0.57827 0.45715
## Proportion of Variance 0.9997 0.00009 0.00005 0.00004 0.00004 0.00003 0.00002
## Cumulative Proportion 0.9997 0.99981 0.99986 0.99989 0.99993 0.99996 0.99998
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.40376 0.1717 2.268e-15 1.854e-15 1.003e-15 8.038e-16
## Proportion of Variance 0.00002 0.0000 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion 1.00000 1.0000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
## PC14
## Standard deviation 4.439e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00Al hacer un summary de las variables vemos que la proporción de varianza explicada por la primera componente es de un 99%, es decir, que es la única relevante. Si nos fijamos vemos que la PC1 está explicada casi al 100% por los m2. En nuestra tabla vemos que los m2 tienen valores mucho más altos que el resto de variables. Por lo tanto, lo que hacemos a continuación es un análisis pero escalando las variables.
prcomp(dummy_data, scale = T)## Standard deviations (1, .., p=14):
## [1] 1.788828e+00 1.554767e+00 1.391248e+00 1.320487e+00 1.255039e+00
## [6] 1.218340e+00 9.006738e-01 7.679348e-01 4.930742e-01 2.297198e-14
## [11] 4.169026e-15 2.788355e-15 1.325148e-15 8.186691e-16
##
## Rotation (n x k) = (14 x 14):
## PC1 PC2 PC3 PC4
## rooms -0.32118124 0.185174155 0.19374853 -0.18231760
## m2 -0.39188093 0.178866531 0.28992752 -0.15926158
## exterior1 -0.22525556 0.323113180 -0.41045579 0.16000993
## exterior0 0.22525556 -0.323113180 0.41045579 -0.16000993
## garage1 -0.39123327 0.104563202 0.02079773 0.20853492
## garage0 0.39123327 -0.104563202 -0.02079773 -0.20853492
## elevator1 -0.28418360 -0.450986811 -0.07378472 0.27110372
## elevator0 0.28418360 0.450986811 0.07378472 -0.27110372
## precio_distritobarato 0.20341246 0.208545636 -0.33198104 0.01261772
## precio_distritomedio -0.02774634 -0.002015776 0.03875669 0.15703618
## precio_distritocaro -0.18188148 -0.213762808 0.30358073 -0.17520956
## tipo_casacasa -0.14205765 0.358486679 0.29646454 -0.20249803
## tipo_casapiso -0.16858177 -0.251362769 -0.39199358 -0.48519863
## tipo_casasotano 0.22431380 0.124024528 0.29030539 0.56945631
## PC5 PC6 PC7 PC8
## rooms 0.17534372 -0.091653082 0.472849806 0.496683083
## m2 0.12920863 -0.038990789 0.209138336 0.106709918
## exterior1 0.31433452 -0.162497819 -0.145528928 -0.064746601
## exterior0 -0.31433452 0.162497819 0.145528928 0.064746601
## garage1 -0.40131726 0.276379019 -0.203710793 0.095975695
## garage0 0.40131726 -0.276379019 0.203710793 -0.095975695
## elevator1 0.11627307 0.012675714 0.278402861 -0.212854201
## elevator0 -0.11627307 -0.012675714 -0.278402861 0.212854201
## precio_distritobarato -0.04064647 0.512859981 0.434958282 -0.014241486
## precio_distritomedio -0.38476634 -0.695032433 0.035559829 0.020880654
## precio_distritocaro 0.43936573 0.186858173 -0.486898651 -0.006820754
## tipo_casacasa -0.08089606 -0.001818982 0.145439408 -0.763400011
## tipo_casapiso -0.15088701 -0.039033227 -0.059791130 0.108495318
## tipo_casasotano 0.18370301 0.040471673 0.007300261 0.171061449
## PC9 PC10 PC11 PC12
## rooms -0.53136448 1.871326e-14 6.573067e-16 -8.713797e-16
## m2 0.79477851 5.562477e-15 5.135097e-16 5.962879e-16
## exterior1 -0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## exterior0 0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## garage1 -0.06271669 -7.002484e-04 -1.864879e-03 2.225429e-02
## garage0 0.06271669 -7.002484e-04 -1.864879e-03 2.225429e-02
## elevator1 -0.02034701 -1.325824e-03 -1.363473e-02 3.145845e-02
## elevator0 0.02034701 -1.325824e-03 -1.363473e-02 3.145845e-02
## precio_distritobarato 0.09621994 -7.055211e-02 4.768696e-02 -5.771886e-01
## precio_distritomedio 0.02039646 -7.038638e-02 4.757494e-02 -5.758327e-01
## precio_distritocaro -0.12064810 -6.816645e-02 4.607446e-02 -5.576714e-01
## tipo_casacasa -0.22215925 -2.483203e-01 -1.912383e-03 3.007261e-02
## tipo_casapiso 0.04299775 -6.861954e-01 -5.284580e-03 8.310109e-02
## tipo_casasotano 0.03813232 -6.729672e-01 -5.182706e-03 8.149911e-02
## PC13 PC14
## rooms -4.757936e-17 3.988302e-17
## m2 -1.440589e-18 -3.087024e-16
## exterior1 1.099278e-02 3.536496e-04
## exterior0 1.099278e-02 3.536496e-04
## garage1 2.157713e-02 -7.064242e-01
## garage0 2.157713e-02 -7.064242e-01
## elevator1 -7.059752e-01 -2.053508e-02
## elevator0 -7.059752e-01 -2.053508e-02
## precio_distritobarato -2.595459e-02 -1.903173e-02
## precio_distritomedio -2.589362e-02 -1.898702e-02
## precio_distritocaro -2.507696e-02 -1.838819e-02
## tipo_casacasa 1.806856e-03 1.253756e-03
## tipo_casapiso 4.992972e-03 3.464565e-03
## tipo_casasotano 4.896719e-03 3.397776e-03Tras escalar las variables podemos observar que están mucho más repartidos los pesos en las componentes principales debido a que estás asegurando que todas las variables tengan la misma importancia en el análisis, independientemente de sus escalas originales y unidades de medida. Con esto evitas el dominio de variables con mayor magnitud y mejoras la interpretación de resultados.
pca_result <- prcomp(dummy_data, center = TRUE, scale = TRUE)
summary(pca_result)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.7888 1.5548 1.3912 1.3205 1.2550 1.2183 0.90067
## Proportion of Variance 0.2286 0.1727 0.1383 0.1245 0.1125 0.1060 0.05794
## Cumulative Proportion 0.2286 0.4012 0.5395 0.6640 0.7765 0.8826 0.94051
## PC8 PC9 PC10 PC11 PC12 PC13
## Standard deviation 0.76793 0.49307 2.297e-14 4.169e-15 2.788e-15 1.325e-15
## Proportion of Variance 0.04212 0.01737 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion 0.98263 1.00000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
## PC14
## Standard deviation 8.187e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00Observamos que con PC6 ya llegamos al 88% de la proporción de varianza explicada pero vemos algunos saltos interesantes. Para ver cuantas componentes explican el modelo hacemos la regla de codo, método heurístico para determinar el número óptimo de clusters en un conjunto de datos.
Al observar el gráfico vemos que la forma del codo sería en la 6|7 aunque el gráfico no es muy claro.
summary_pca <- summary(pca_result)
var_exp <- summary_pca$importance[2,]
barplot(var_exp, col = "#0BB363",
main = "Gráfico PCA - Variance explained",
xlab = "Principal Components",
ylab = "Proportion of Variance Explained")
lines(var_exp, col = "#DC143C", type = "l", lwd = 2)
points(var_exp, col = "#DC143C", pch = 20, cex = 1.2)
A modo de ejercicio para comprender el funcionamiento de PCA puede resultar útil hacer estos cálculos, sin embargo para el dataset en cuestión ya se pueden formar modelos robustos con muy pocas variables (los m2 de superficie, por ejemplo, ya son suficientes para predecir un gran porcentaje de precios).
10 Aprendizaje supervisado
10.1 GLM
En aprendizaje automático, GLM (Generalized Linear Models) se refiere a una clase de modelos lineales que se utilizan para predecir una variable respuesta a partir de una o más variables explicativas. Tiene tres componentes: Variable respuesta, función de enlace y función lineal sistemática, y son muy útiles para utilizarse en clasificación, regresión y problemas de predicción.
Aplicado a nuestro dataset, lo utilizaremos para predecir el precio binario de las casas (cara/barata).
train_glm <- mhp_train[,c(3:7,10,11,12)]
test_glm <- mhp_test[,c(3:7,10,11,12)]train_glm[,2:3]<- scale(train_glm[,2:3])
test_glm[,2:3]<- scale(test_glm[,2:3])Tras seleccionar las variables y normalizarlas para que todas tengan el mismo peso, desarrollamos el modelo con la variable de respuesta “precio_bin” y las variables independientes “exterior”, “garage”, “elevator”, “rooms”, “precio_distrito” y “tipo_casa”, usando una distribución binomial.
glm_mhp = glm(formula = precio_bin ~ exterior + garage + elevator + rooms + m2 + precio_distrito + tipo_casa,
data = train_glm,
family = binomial)
summary(glm_mhp)##
## Call:
## glm(formula = precio_bin ~ exterior + garage + elevator + rooms +
## m2 + precio_distrito + tipo_casa, family = binomial, data = train_glm)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5396 -0.1705 0.0219 0.2836 5.0687
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.99261 0.62333 -1.592 0.11129
## exterior0 -0.35269 0.13261 -2.660 0.00782 **
## garage0 0.79113 0.10737 7.368 1.73e-13 ***
## elevator0 1.67266 0.12813 13.054 < 2e-16 ***
## rooms -0.01548 0.08151 -0.190 0.84940
## m2 -6.20689 0.23649 -26.246 < 2e-16 ***
## precio_distritomedio -2.52394 0.11266 -22.403 < 2e-16 ***
## precio_distritocaro -4.33905 0.14476 -29.974 < 2e-16 ***
## tipo_casapiso 0.34751 0.62328 0.558 0.57716
## tipo_casasotano 1.01944 0.62491 1.631 0.10282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10718.8 on 7731 degrees of freedom
## Residual deviance: 3535.9 on 7722 degrees of freedom
## AIC: 3555.9
##
## Number of Fisher Scoring iterations: 8Tras realizar muchas comprobaciones y probar de maneras diversa el modelo no conseguimos interpretar bien los datos. No somos capaces de entender por qué sale inversamente proporcional los m2,el precio distrito mdio y el caro, que precisamente son las 3 variables que más definen que una casa sea cara, entonces no entendemos qué significa en este caso que salgan negativas. Lo dejamos pendiente para que por favor nos expliqueis a qué se puede deber esto y así saberlo para futuras ocasiones.
head(predict(glm_mhp))## 1 2 3 4 5 6
## -2.6264925 1.5662145 -6.1858748 0.1476553 0.4424557 -3.7636044Esto sirve para predecir la probabilidad logit ajustado previamente y para mostrar las primeras seis predicciones. Cada número representa la probabilidad logit de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.
ajustados <- fitted(glm_mhp)
head(fitted(glm_mhp))## 1 2 3 4 5 6
## 0.067452746 0.827243283 0.002054072 0.536846908 0.608844012 0.022673932Aquí obtenemos las probabilidades ajustadas que representan las probabilidades estimadas de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.
Gráficos de residuales: Estos gráficos permiten evaluar la calidad del ajuste del modelo. Puedes trazar residuales de Pearson, deviance, o residuales estandarizados frente a las predicciones ajustadas o valores ajustados para examinar si hay patrones en los residuales. Un buen ajuste del modelo no mostrará patrones claros en los residuales.
residuales <- residuals(glm_mhp, type = "pearson")
head(residuales)## 1 2 3 4 5 6
## -0.26894557 -2.18826117 -0.04536849 -1.07662124 0.80153404 -0.15231535data_residuales <- data.frame(Residuales = residuales, Ajustados = ajustados)
ggplot(data_residuales, aes(x = Ajustados, y = Residuales)) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(title = "Gráfico de residuales", x = "Valores ajustados", y = "Residuales") + ylim(-10,10)
predicciones_glm <- ifelse(test = glm_mhp$fitted.values > 0.5, no = 0, yes = 1)
matriz_confusion_glm <- table(glm_mhp$model$precio_bin, predicciones_glm,
dnn = c("observaciones", "predicciones"))
matriz_confusion_glm## predicciones
## observaciones 0 1
## 1 3449 414
## 0 329 3540Finalmente generamos predicciones binarias basadas en las probabilidades ajustadas del modelo lineal generalizado y creamos una matriz de confusión para evaluar el rendimiento del modelo donde se muestra la distribución de las predicciones correctas e incorrectas en función de las observaciones reales y las predichas.
Obersvamos que el modelo no acierta nada, esto imagino que estará relacionado con lo anterior que no llegamos a entender.
La matriz se lee de la siguiente manera:
Verdaderos negativos (VN): La esquina superior izquierda representa el número de observaciones clasificadas correctamente como “0” (clase negativa).
Falsos positivos (FP): La esquina superior derecha representa el número de observaciones que en realidad eran “0” pero fueron clasificadas incorrectamente como “1” (clase positiva) por el modelo.
Falsos negativos (FN): La esquina inferior izquierda representa el número de observaciones que en realidad eran “1” pero fueron clasificadas incorrectamente como “0” por el modelo.
Verdaderos positivos (VP): La esquina inferior derecha representa el número de observaciones clasificadas correctamente como “1”.
10.2 KNN
El modelo KNN (K-Nearest Neighbors) es un algoritmo de aprendizaje supervisado que se utiliza tanto para clasificación como para regresión. Es un método basado en instancias que no hace supuestos sobre la forma funcional de la relación entre las variables independientes y la variable dependiente. KNN utiliza la información de los vecinos más cercanos para realizar predicciones. Para ello calcula distancias como la Euclidiana, Manhattan, Minkowski, etc…
La idea detrás de KNN es simple, una observación se clasifica o se le asigna un valor basado en las características de sus K vecinos más cercanos.
Algunas ventajas del algoritmo KNN son su simplicidad o su capacidad para manejar relaciones no lineales entre variables. Sin embargo también tiene algunas desventajas, como su sensibilidad a la maldición de la dimensionalidad, el efecto de las variables irrelevantes y el costo computacional asociado con el cálculo de distancias en conjuntos de datos grandes.
Procedemos al cálculo de este modelo e intentamos ver si hay un número correcto de k-vecinos que podemos usar.
long = 15
accuracy = rep(0,long)
f1score = rep(0,long)
recall = rep(0,long)
precision = rep(0,long)
for (i in 1:long)
{
prediccion_knn_cv =knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")],
k=i, cl=mhp_train$precio_bin)
accuracy[i] = sum(prediccion_knn_cv == mhp_train$precio_bin) /nrow(mhp_train)
recall[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & mhp_train$precio_bin == TRUE) / sum(mhp_train$precio_bin == TRUE)
precision[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & prediccion_knn_cv == TRUE) / sum(prediccion_knn_cv == TRUE)
f1score[i] = 2*precision[i]*recall[i]/(precision[i]+recall[i])
}
resultados_knn = as.data.frame(cbind(accuracy,f1score,precision,recall))
resultados_knn = resultados_knn %>% mutate(index=as.factor(seq(1:long)))
max(resultados_knn$f1score)## [1] NaNwhich.max(resultados_knn$f1score)## integer(0)ggplot(data=resultados_knn,aes(x=index,y=accuracy)) +
geom_col(colour="cyan4",fill="cyan3")+
ggtitle("Accuracy")
ggplot(data=resultados_knn,aes(x=index,y=f1score)) +
geom_col(colour="orange4",fill="orange3") +
ggtitle("F1_score values")
Tratamos de encontrar el valor óptimo de k utilizando validación cruzada. Se calculan varias métricas de evaluación, como precisión, recall y F1-score, para cada valor de k y finalmente tratamos de graficar la precisión en función del valor de k, pero analizando el gráfico no observamos ninguna diferencia, así que cogemos 5, que es el tamaño por convención que se suele coger.
Pasamos a realizar predicciones en el conjunto de datos de entrenamiento y en el de test además de realizar la matriz de confusión.
# En train
prediccion_knn5_train <- knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")],
k=5, cl=mhp_train$precio_bin, prob = TRUE)
confusionMatrix(table(prediccion_knn5_train,mhp_train$precio_bin), positive= "1")## Confusion Matrix and Statistics
##
##
## prediccion_knn5_train 1 0
## 1 3233 531
## 0 630 3338
##
## Accuracy : 0.8498
## 95% CI : (0.8417, 0.8577)
## No Information Rate : 0.5004
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6997
##
## Mcnemar's Test P-Value : 0.004026
##
## Sensitivity : 0.8369
## Specificity : 0.8628
## Pos Pred Value : 0.8589
## Neg Pred Value : 0.8412
## Prevalence : 0.4996
## Detection Rate : 0.4181
## Detection Prevalence : 0.4868
## Balanced Accuracy : 0.8498
##
## 'Positive' Class : 1
## #En test
prediccion_knn5_test <- knn(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], mhp_test[,c("exterior","rooms","m2","elevator", "garage")],
k=5, cl=mhp_train$precio_bin, prob = TRUE)
confusionMatrix(table(prediccion_knn5_test,mhp_test$precio_bin), positive= "1")## Confusion Matrix and Statistics
##
##
## prediccion_knn5_test 1 0
## 1 1577 250
## 0 356 1683
##
## Accuracy : 0.8432
## 95% CI : (0.8314, 0.8546)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.6865
##
## Mcnemar's Test P-Value : 1.996e-05
##
## Sensitivity : 0.8158
## Specificity : 0.8707
## Pos Pred Value : 0.8632
## Neg Pred Value : 0.8254
## Prevalence : 0.5000
## Detection Rate : 0.4079
## Detection Prevalence : 0.4726
## Balanced Accuracy : 0.8432
##
## 'Positive' Class : 1
## La matriz de confusión en train nos muestra la cantidad de predicciones correctas e incorrectas de cada clase (1 y 0) en el conjunto de entrenamiento.
Accuracy (exactitud): 0.8482. Es la proporción de predicciones correctas en el conjunto de entrenamiento. En otras palabras, el modelo KNN clasifica correctamente el 84.82% de las observaciones del conjunto de entrenamiento.
Kappa: 0.6963. El índice kappa mide la concordancia entre las predicciones y los valores reales. Un kappa de 1 indica una concordancia perfecta, mientras que un kappa de 0 indica una concordancia totalmente aleatoria. Nuestro kappa de 0.6963 muestra una buena concordancia entre las predicciones y los valores reales en el conjunto de entrenamiento.
Mcnemar’s Test P-Value: 0.003201. Esta prueba evalúa si hay una diferencia significativa entre las predicciones falsas positivas y las predicciones falsas negativas.
La matriz de confusión en test nos muestra valores muy parecidos en los indicadores anteriores lo que nos lleva a pensar que el modelo tiene un rendimiento consistente y generaliza bien lo datos ya que el modelo clasifica correctamente las observaciones y no parece sobreajustado.
Hemos realizado diferentes pruebas quitando variables y se observa que si quitamos la variable m2 nos dice que hay muchos empates. Sin embargo, quitando las otras variables el acierto en train y test varía.
Obteniendo más o menos un 84% de acierto creemos que nos puede servir para clasificar las nuevas casas que entrasen en el dataset, puesto que no hemos sido capaces de mejorar ese %.
10.3 Decision Trees
El modelo de Árbol de Decisión es un algoritmo de aprendizaje supervisado que funciona dividiendo iterativamente el conjunto de datos en subconjuntos basados en las características del mismo. Al final de este proceso, se genera un árbol con nodos y hojas que representan las divisiones y las decisiones.
Algunas características importantes de los árboles de decisión son la interpretabilidad, capacidad para manejar datos no lineales, el poco preprocesamiento de los datos, la sensibilidad a datos desequilibrados y el posible sobreajuste en el caso de hacer un arbol con demasiadas ramificaciones.
A continuación hacemos el modelo de arbol de decisión.
arbol <- rpart(precio_bin ~ garage + elevator + rooms + tipo_casa + precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")
Para evaluar el modelo utilizamos la función predict sobre el conjunto test usando el argumento type = “class” que indica que se espera una salida categórica creandose una tabla de contingencias “tab1” con las predicciones que coinciden o no respecto a las observaciones reales.
Calculamos y obtenemos el número 0.8396275 que representa la tasa de acierto (accuracy) del modelo. El Árbol de Decisión hizo predicciones correctas en casi el 84% de los casos en el conjunto de datos de prueba. Nos parece un porcentaje bastante aceptable.
tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1## [1] 0.8396275Siguiendo con la validación del modelo podemos usar las funciones printcp y plotcp que se utilizan para analizar la complejidad y el rendimiento usando el “pruning” o poda de ramas del arbol que a menudo reduce su complejidad y permite una mejor generalización del mismo disminuyendo el subreajuste.
La función printcp es una función que muestra una tabla que contiene información sobre la complejidad y la tasa de error en diferentes niveles de tamaño del árbol proporcionando una visión general de cómo cambia el rendimiento a medida que aumenta la complejidad. El valor de cp ha de ser tal que la tasa de error de validación cruzada sea mínima.
La función plotcp crea un gráfico de la tasa de error. Este gráfico es útil para visualizar cómo se comporta el rendimiento del árbol a medida que aumenta su tamaño y ayuda a decidir qué tamaño de árbol es el más apropiado para el problema en cuestión.
printcp(arbol)##
## Classification tree:
## rpart(formula = precio_bin ~ garage + elevator + rooms + tipo_casa +
## precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
##
## Variables actually used in tree construction:
## [1] elevator precio_distrito rooms
##
## Root node error: 3863/7732 = 0.49961
##
## n= 7732
##
## CP nsplit rel error xerror xstd
## 1 0.516438 0 1.00000 1.03055 0.0113762
## 2 0.044784 1 0.48356 0.48356 0.0097435
## 3 0.041677 3 0.39399 0.40098 0.0091108
## 4 0.016179 4 0.35232 0.35232 0.0086689
## 5 0.010000 6 0.31996 0.31996 0.0083418Obtenidos los resultados, pasamos a analizarlos:
- CP (Cost Complexity): Mide la complejidad asociado con cada nivel de poda en el árbol de decisión. Un CP más pequeño indica una complejidad menor.
- nsplit: Número de divisiones (splits) realizadas hasta el momento en el árbol.
- rel error: Error relativo en cada etapa de poda del árbol.
- xerror: Error de validación cruzada (cross-validated) estimado para el árbol en cada etapa de poda.
- xstd: Desviación estándar del error de validación cruzada (xerror).
En base a los resultados presentados, se puede observar que el árbol sin poda tiene un error relativo de 1.0 y un error de validación cruzada de 1.0189. A medida que se realizan más divisiones en el árbol el CP disminuye y también lo hace el error relativo y el error de validación cruzada.
Para seleccionar el árbol óptimo, generalmente se busca el CP más pequeño que tenga un error de validación cruzada dentro de una desviación estándar del mínimo error de validación cruzada.
arbol$cptable[which.min(arbol$cptable[, "xerror"]),
"CP"]## [1] 0.01El propósito es identificar el mejor árbol de decisión con respecto a su capacidad de generalización en datos no vistos. La función which.min() devuelve el índice del valor mínimo en la columna “xerror” de la tabla de complejidad cptable del árbol de decisión. Luego, con ese índice, se extrae el valor correspondiente de la columna “CP” en la tabla de complejidad.
plotcp(arbol)
Podamos el árbol como indica la función.
pruneTREE1 = prune(arbol, cp = arbol$cptable[which.min(arbol$cptable[,
"xerror"]), "CP"])
fancyRpartPlot(pruneTREE1, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")
pruneTREE2 = prune(arbol, cp = arbol$cptable[5, "CP"])
fancyRpartPlot(pruneTREE2, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")
tab2 = table(pred = predict(pruneTREE2, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
acierto2 = sum(diag(tab2))/ntest
acierto2## [1] 0.8396275Tras la poda nos sale literalmente es mismo resultado.
Creamos la matriz de confusión.
prediccion_1 <- predict(arbol, newdata = mhp_train, type = "class")confusionMatrix(prediccion_1, mhp_train[["precio_bin"]])## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 0
## 1 3246 619
## 0 617 3250
##
## Accuracy : 0.8401
## 95% CI : (0.8318, 0.8482)
## No Information Rate : 0.5004
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.6803
##
## Mcnemar's Test P-Value : 0.9773
##
## Sensitivity : 0.8403
## Specificity : 0.8400
## Pos Pred Value : 0.8398
## Neg Pred Value : 0.8404
## Prevalence : 0.4996
## Detection Rate : 0.4198
## Detection Prevalence : 0.4999
## Balanced Accuracy : 0.8401
##
## 'Positive' Class : 1
## Los resultados en la matriz de confusión son los siguientes:
- Verdaderos positivos (VP): 3246
- Verdaderos negativos (VN): 3250
- Falsos positivos (FP): 619
- Falsos negativos (FN): 617
También obtenemos:
- Accuracy (exactitud): 84.01% de predicciones correctas.
- Kappa: 68.03%
- Sensitivity (sensibilidad): 0.8403. También conocida como recall o tasa verdadera positiva, mide la proporción de casos positivos que se identificaron correctamente.
- Specificity (especificidad): 0.8400. Mide la proporción de casos negativos que se identificaron correctamente. Un valor más alto es mejor.
- Pos Pred Value (valor predictivo positivo): 0.8398. Mide la proporción de predicciones positivas que son realmente positivas. Un valor más alto es mejor.
- Neg Pred Value (valor predictivo negativo): 0.8404. Mide la proporción de predicciones negativas que son realmente negativas. Un valor más alto es mejor.
- Mcnemar’s Test P-Value: 0.9773. Esta prueba compara las diferencias entre los falsos positivos y los falsos negativos. Un valor alto de “p” (como en este caso) sugiere que no hay una diferencia significativa entre ellos.
En resumen, este modelo de árbol de decisión podado tiene una exactitud del 84.01% y un coeficiente kappa de 0.6803, lo que indica un buen rendimiento en la clasificación.
Sin embargo para este modelo se han tenido en cuenta bastantes variables, y aunque el rendimiento es bueno, se puede alcanzar un rendimiento casi igual utilizando sólo los metros cuadrados, como se puede ver a continuación
arbol <- rpart(precio_bin ~ m2, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")
tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1## [1] 0.8344542Esta es una forma simple de confirmar la relevancia que tiene la superficie de una casa sobre su precio final, en contraste con otras variables como el tener ascensor o garaje. Además, esto nos deja un modelo mucho más simple y explicable.
10.4 Random Forest
Algoritmo de aprendizaje supervisado que se basa en la construcción de múltiples árboles de decisión durante el entrenamiento y la combinación de sus resultados para hacer una predicción. El objetivo principal de este enfoque es mejorar la precisión y la robustez del modelo en comparación con un único árbol de decisión.
El funcionamiento de un algoritmo Random Forest sería el siguiente:
Seleccionar muestras de arranque: El algoritmo selecciona aleatoriamente un subconjunto de muestras del conjunto de datos de entrenamiento, con reemplazo. Esto significa que una muestra puede ser seleccionada más de una vez.
Crear un árbol de decisión: Se crea un árbol de decisión utilizando el subconjunto de muestras seleccionado. Durante la construcción del árbol, en cada división del nodo, se selecciona un subconjunto aleatorio de características como candidatas para la división. La mejor división se elige en función de la medida de impureza como la entropía o el índice de Gini. Esto introduce aleatoriedad adicional en el proceso, lo que ayuda a reducir la correlación entre los árboles.
Repetir: El proceso se repite varias veces (número de árboles especificado) para construir un conjunto de árboles de decisión.
Combinar las predicciones: Para hacer una predicción en un nuevo caso, el algoritmo Random Forest ejecuta todos los árboles individualmente y luego combina sus resultados. En el caso de la clasificación, la clase con más votos (predicciones) de todos los árboles se elige como la predicción final.
Pasamos a construir el modelo comenzando con 10 árboles aleatórios.
# Escalamos la columna 4 en train y test
set.seed(19042023)
mhp_train_scaled <- mhp_train
mhp_train_scaled[, 4] <- scale(mhp_train[, 4])
mhp_test_scaled <- mhp_test
mhp_test_scaled[, 4] <- scale(mhp_test[, 4])
# Nos quedamos con todas las variables menos barrio y distrito porque hay muchas categorias
# Cogemos solo 10 árboles en el parámetro ntree
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
y = mhp_train_scaled$precio_bin,
ntree = 10)# Predicción de los resultados con el conjunto de testing
y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 0
## 1 1676 127
## 0 257 1806
##
## Accuracy : 0.9007
## 95% CI : (0.8908, 0.9099)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.8013
##
## Mcnemar's Test P-Value : 4.61e-11
##
## Sensitivity : 0.8670
## Specificity : 0.9343
## Pos Pred Value : 0.9296
## Neg Pred Value : 0.8754
## Prevalence : 0.5000
## Detection Rate : 0.4335
## Detection Prevalence : 0.4664
## Balanced Accuracy : 0.9007
##
## 'Positive' Class : 1
## Ahora probamos con 100 árboles aleatorios para ver si cambia mucho el resultado.
set.seed(19042022)
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
y = mhp_train_scaled$precio_bin,
ntree = 100)y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])## Confusion Matrix and Statistics
##
## Reference
## Prediction 1 0
## 1 1684 108
## 0 249 1825
##
## Accuracy : 0.9077
## 95% CI : (0.8981, 0.9166)
## No Information Rate : 0.5
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.8153
##
## Mcnemar's Test P-Value : 1.267e-13
##
## Sensitivity : 0.8712
## Specificity : 0.9441
## Pos Pred Value : 0.9397
## Neg Pred Value : 0.8799
## Prevalence : 0.5000
## Detection Rate : 0.4356
## Detection Prevalence : 0.4635
## Balanced Accuracy : 0.9077
##
## 'Positive' Class : 1
## La precisión del modelo es del 0,9077 y el intervalo de confianza del 95% es (0,8908, 0,9099), lo que indica que podemos estar bastante seguros de que la precisión real del modelo está dentro de este rango. El valor de kappa de 0,8013 indica que hay una buena concordancia entre las predicciones del modelo y las observaciones reales.
Además, la sensibilidad del modelo es del 0,8670 y la especificidad del 0,9343. La tasa de detección es del 0,4335 y la prevalencia es del 0,5000. En general, el modelo tiene un desempeño bastante bueno en la clasificación de las viviendas por encima y por debajo del umbral.
10.5 SVM
Support Vector Machine, algoritmo de aprendizaje supervisado que se utiliza tanto para problemas de clasificación como de regresión. En el caso de la clasificación, el objetivo principal de SVM es encontrar el hiperplano óptimo que separa los datos en diferentes clases. Un hiperplano es un límite de decisión que divide el espacio de características en dos partes, de modo que cada parte contenga puntos de datos de una clase específica. En el caso de la regresión, SVM busca encontrar una función que tenga el menor error de ajuste posible.
El algoritmo SVM funciona de la siguiente manera:
Encuentra el hiperplano óptimo: El algoritmo busca el hiperplano que tenga el margen máximo entre las dos clases. El margen es la distancia entre el hiperplano y los puntos de datos más cercanos a él, llamados “vectores de soporte”. Estos vectores de soporte son los puntos de datos que definen el límite de decisión y son fundamentales para el proceso de SVM.
Solución de un problema de optimización: La búsqueda del hiperplano óptimo implica la solución de un problema de optimización convexa. Este problema se resuelve mediante técnicas de optimización como el método de Lagrange o el algoritmo de optimización secuencial mínima (SMO).
Transformación a través del kernel: En casos en los que los datos no son linealmente separables en el espacio de características original, SVM puede utilizar una función de kernel para transformar los datos en un espacio de características de mayor dimensión donde se puedan separar linealmente. Los kernels comunes incluyen el kernel lineal, el kernel polinómico, el kernel de base radial (RBF) y el kernel sigmoide. La elección del kernel adecuado es crucial para el rendimiento del modelo SVM.
Clasificación: Para hacer una predicción en un nuevo punto de datos, el algoritmo SVM evalúa la posición del punto con respecto al hiperplano óptimo. Si el punto está por encima del hiperplano, se asigna a una clase, y si está por debajo, se asigna a la otra clase.
Para resolver SVM se necesita conocer el kernel y el parámetro de regularización, que es una cuota superior sobre las variables duales, en definitiva, un coste. Pasamos a construirun primer el modelo usando un kernel radial RBF.
# Estandarizamos las variables numéricas "m2" y "rooms"
mhp_train_svm <- mhp_train
mhp_train_svm$m2 <- scale(mhp_train$m2)
mhp_train_svm$rooms <- scale(mhp_train$rooms)
mhp_test_svm <- mhp_test
mhp_test_svm$m2 <- scale(mhp_test$m2)
mhp_test_svm$rooms <- scale(mhp_test$rooms)
# Entrenamos el modelo SVM kernel con kernel radial
modelo_svm <- svm(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito, data = mhp_train_svm, kernel = "radial", cost = 10, scale = TRUE)
# Hacemos predicciones en el conjunto de prueba
svm_predicciones <- predict(modelo_svm, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")])
# Calculamos la matriz de confusión y la precisión del modelo
svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test_svm$precio_bin)
print(svm_matriz_confusion)## Real
## Prediccion 1 0
## 1 1692 116
## 0 241 1817precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
print(paste0("Precisión del modelo SVM: ", round(precision * 100, 2), "%"))## [1] "Precisión del modelo SVM: 90.77%"Esto significa que:
- 1692 casas con precio alto fueron clasificadas correctamente como precio alto (verdaderos positivos).
- 1817 casas con precio bajo fueron clasificadas correctamente como precio bajo (verdaderos negativos).
- 241 casas con precio alto fueron clasificadas incorrectamente como precio bajo (falsos negativos).
- 116 casas con precio bajo fueron clasificadas incorrectamente como precio alto (falsos positivos).
La precisión del modelo es del 90.77%.
Elección del Kernel
La elección del kernel en un modelo SVM depende de la naturaleza y la distribución de los datos. No hay una regla única para seleccionar el kernel, pero aquí hay algunas pautas generales para ayudarte a tomar una decisión:
Kernel lineal: Si tus datos son linealmente separables o la cantidad de atributos es alta en comparación con la cantidad de muestras, entonces un kernel lineal suele funcionar bien. Además, es computacionalmente más eficiente que otros kernels.
Kernel polinomial: Si tus datos tienen una estructura que puede ser capturada por un polinomio de un grado específico, el kernel polinomial podría ser una buena opción. Sin embargo, a medida que aumenta el grado del polinomio, el riesgo de sobreajuste aumenta y la complejidad computacional puede crecer.
Kernel radial (RBF): El kernel radial es una opción popular y versátil en muchos casos, especialmente cuando no se sabe de antemano si los datos son linealmente separables. Puede manejar datos no lineales y complejos. Sin embargo, encontrar los parámetros adecuados (costo y gamma) puede requerir una búsqueda en el espacio de parámetros y, en consecuencia, un mayor tiempo de cómputo.
Kernel sigmoide: El kernel sigmoide es similar a una función de activación de una red neuronal y puede ser una opción si deseas un enfoque similar al de una red neuronal pero con un modelo SVM.
Vamos a probar a representar los datos para ver si de alguna manera puede tomarse una mejor decisión sobre el Kernel a usar.
mhp_train %>%
ggplot(aes(x = m2, y = rooms, color = as.factor(precio_bin))) +
geom_point() +
facet_grid(elevator ~ garage, labeller = labeller(
elevator = c('0' = 'Sin Elevador', '1' = 'Con Elevador'),
garage = c('0' = 'Sin Garaje', '1' = 'Con Garaje')
)) +
labs(title = "Diagrama de dispersión de m2 vs rooms, según garage y elevator",
x = "Metros cuadrados",
y = "Número de habitaciones",
color = "Precio bin") +
scale_color_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
plot.title.position = "plot")
Probamos varios Kernels.
Creamos una función entrenar_evaluar_svm() que entrenará y evaluará un modelo SVM utilizando el tipo de kernel especificado en el argumento kernel_type. Luego con lapply() se aplicará esta función a cada tipo de kernel en el vector kernels. La función presentará un matriz de confusión y la precisión para cada tipo de kernel.
# Definir una función para entrenar y evaluar un modelo SVM con un tipo de kernel específico
entrenar_evaluar_svm <- function(kernel_type) {
# Entrenamos el modelo SVM con el kernel especificado
modelo_svm <- svm(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito, data = mhp_train_svm, kernel = kernel_type, cost = 10, scale = TRUE)
# Hacemos predicciones en el conjunto de prueba
svm_predicciones <- predict(modelo_svm, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")])
# Calculamos la matriz de confusión y la precisión del modelo
svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test_svm$precio_bin)
precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
# Devolvemos una lista con la precisión y la matriz de confusión
return(list(precision = precision, matriz_confusion = svm_matriz_confusion))
}
# Probar diferentes tipos de kernel
kernels <- c("linear", "polynomial", "radial", "sigmoid")
resultados <- lapply(kernels, entrenar_evaluar_svm)
# Imprimir los resultados
for (i in 1:length(kernels)) {
cat(paste("Precisión del modelo SVM con kernel", kernels[i], ":", round(resultados[[i]]$precision * 100, 2), "%\n"))
print(resultados[[i]]$matriz_confusion)
}## Precisión del modelo SVM con kernel linear : 90.46 %
## Real
## Prediccion 1 0
## 1 1706 142
## 0 227 1791
## Precisión del modelo SVM con kernel polynomial : 90.43 %
## Real
## Prediccion 1 0
## 1 1688 125
## 0 245 1808
## Precisión del modelo SVM con kernel radial : 90.77 %
## Real
## Prediccion 1 0
## 1 1692 116
## 0 241 1817
## Precisión del modelo SVM con kernel sigmoid : 85 %
## Real
## Prediccion 1 0
## 1 1633 280
## 0 300 1653Se observa que el KERNEL RADIAL es el que mejor predice, ahora vamos a ajustar algunos parámetros del modelo utilizando la función tune para encontrar los mejores parámetros.
# Entrenar un modelo SVM con diferentes parámetros
modelo_svm_radial <- svm(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito, data = mhp_train_svm, kernel = "radial", cost = 10, scale = TRUE)
# Utilizar la función tune() para encontrar los mejores parámetros
tuned_parameters <- tune(svm,
precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito,
data = mhp_train_svm,
kernel = "radial",
ranges = list(cost = c(0.1, 1, 10, 100, 1000),
gamma = c(0.1, 0.5, 1, 2, 5)))
# Entrenar un modelo SVM con los mejores parámetros encontrados
best_model <- tuned_parameters$best.model
# Hacer predicciones en el conjunto de prueba con el mejor modelo
best_predictions <- predict(best_model, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")])
# Calcular la matriz de confusión y la precisión del mejor modelo
best_matriz_confusion <- table(Prediccion = best_predictions, Real = mhp_test_svm$precio_bin)
print(best_matriz_confusion)## Real
## Prediccion 1 0
## 1 1692 118
## 0 241 1815best_precision <- sum(diag(best_matriz_confusion)) / sum(best_matriz_confusion)
print(paste0("Precisión del mejor modelo SVM: ", round(best_precision * 100, 2), "%"))## [1] "Precisión del mejor modelo SVM: 90.71%"La precisión del modelo no mejora del 90.77% anterior y el coste computacional en la busqueda del mejor modelo es enorme. Parece que los parámetros del modelo elegidos durante el proceso de afinación no han sido los correctos. La afinación e hiperparámetros implica buscar la mejor combinación de hiperparámetros dentro del rango de búsqueda que se proporciona. Si el rango de búsqueda no contiene la combinación óptima de hiperparámetros, el proceso de afinación podría terminar eligiendo una combinación subóptima que resulte en un rendimiento inferior.
11 Evaluación y comparación de modelos
Evaluamos y comparamos los modelos que mayor porcentaje de predicción nos han proporcionado, Random Forest y SVM radial cost 10 (modelo inicial que nos daba un alto porcentaje de predicción y bajo coste computacional).
11.1 Random Forest
11.1.1 K-fold Cross Validation
K-fold Cross Validation es un método de validación cruzada donde la idea principal es dividir el conjunto de datos en K subconjuntos más pequeños, llamados “folds”, y utilizarlos para entrenar y probar el modelo de manera iterativa. Esto ayuda a evitar el sobreajuste y proporciona una estimación más precisa del rendimiento del modelo en datos no vistos.
El proceso de K-fold Cross Validation es el siguiente:
- Dividir el conjunto de datos en K subconjuntos (folds) de
aproximadamente el mismo tamaño.
- Para cada fold:
- Utilizar K-1 folds para entrenar el modelo.
- Utilizar el fold restante como conjunto de prueba para evaluar el rendimiento del modelo.
- Para cada fold:
- Repetir el proceso K veces, asegurándose de que cada fold se utilice como conjunto de prueba exactamente una vez.
- Calcular métricas de rendimiento, como precisión, error cuadrático medio, etc., para cada una de las K iteraciones y promediar los resultados para obtener una estimación más sólida del rendimiento del modelo.
La ventaja principal de utilizar K-fold Cross Validation es que todos los datos se utilizan tanto para entrenamiento como para pruebas, lo que proporciona una mejor estimación del rendimiento del modelo. Además, al repetir el proceso K veces, se reduce la varianza en las métricas de rendimiento, lo que resulta en una evaluación más confiable del modelo.
# Renombramos los niveles de la variable de respuesta en los conjuntos de datos de entrenamiento y prueba
mhp_train_scaled$precio_bin <- as.factor(mhp_train_scaled$precio_bin)
levels(mhp_train_scaled$precio_bin) <- c("Level0", "Level1")
mhp_test_scaled$precio_bin <- as.factor(mhp_test_scaled$precio_bin)
levels(mhp_test_scaled$precio_bin) <- c("Level0", "Level1")
# Configuramos la validación cruzada K-fold (por ejemplo, K=10)
k_folds <- 10
folds <- createFolds(mhp_train_scaled$precio_bin, k = k_folds)
# Ajustamos un modelo Random Forest usando validación cruzada K-fold con la librería caret
control <- trainControl(method = "cv", number = k_folds, index = folds)
K_F_modelo_rf <- train(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10, 11)],
y = mhp_train_scaled$precio_bin,
method = "rf",
ntree = 10,
trControl = control)
# Imprimimos el resultado del modelo
print(K_F_modelo_rf)## Random Forest
##
## 7732 samples
## 7 predictor
## 2 classes: 'Level0', 'Level1'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 773, 773, 772, 774, 773, 773, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.8911737 0.7823444
## 4 0.8847361 0.7694716
## 7 0.8757834 0.7515664
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.Estos resultados muestran el rendimiento del modelo de Random Forest para diferentes valores del hiperparámetro mtry que es el número de predictores que se utilizan en cada árbol de decisión. Se han probado 3 valores diferentes: 2, 4 y 7. y para cada valor de mtry se muestran precisiones cercanas al 90% como en el modelo inicial. K-fold validation nos indica que el mayor accuracy se obtiene con tan solo dos predictores, no hay necesidad de más.
11.1.2 N-fold Cross Validation
Es otra técnica de validación donde el conjunto de datos se divide en N subconjuntos y se entrena N veces utilizando cada vez un pliegue diferente como conjunto de validación, y los restantes N-1 pliegues como conjunto de entrenamiento.
Al final del proceso, se obtienen N métricas de rendimiento, una por cada iteración. Estas métricas se promedian para obtener una estimación única del rendimiento del modelo. Esta estimación es más fiable que la que se obtiene utilizando únicamente un único conjunto de entrenamiento y validación, ya que reduce el riesgo de que el rendimiento del modelo se vea afectado por una división aleatoria específica del conjunto de datos.
La validación cruzada de N pliegues ayuda a identificar si un modelo está sobreajustado (overfitting) o infraajustado (underfitting).
# Configura el método de control de entrenamiento
control <- trainControl(method = "cv",
number = 5, # N-fold Cross Validation
savePredictions = "final",
classProbs = TRUE)
# Entrena el modelo de Random Forest con N-fold Cross Validation
N_F_modelo_rf <- train(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
y = mhp_train_scaled$precio_bin,
method = "rf",
trControl = control,
tuneGrid = data.frame(mtry = c(2, 4, 7)))
# Muestra los resultados del modelo
print(N_F_modelo_rf)## Random Forest
##
## 7732 samples
## 7 predictor
## 2 classes: 'Level0', 'Level1'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6186, 6185, 6186, 6186, 6185
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.9053296 0.8106515
## 4 0.9046822 0.8093612
## 7 0.8909726 0.7819452
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.Para hacer N-fold Validation se entrenó un modelo de Random Forest utilizando validación cruzada de 5 pliegues, y se encontró que el valor óptimo de mtry para este conjunto de datos y clasificación es 2. El modelo alcanzó una exactitud promedio del 90.4% y un índice Kappa promedio de 0.809.
11.2 SVM
Repetimos las validaciones para SVM.
11.2.1 K-fold Cross Validation
# fórmula
formula_svm <- precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito
# Configura los parámetros de entrenamiento
control_svm <- trainControl(method = "cv", number = 10) # 10-fold Cross Validation
# Entrena el modelo SVM con K-fold Cross Validation
K_F_modelo_svm_cv <- train(formula_svm,
data = mhp_train_svm,
method = "svmRadial",
trControl = control_svm,
preProcess = c("center", "scale"),
tuneGrid = data.frame(sigma = 0.1, C = 10))
# Muestra los resultados
print(K_F_modelo_svm_cv)## Support Vector Machines with Radial Basis Function Kernel
##
## 7732 samples
## 6 predictor
## 2 classes: '1', '0'
##
## Pre-processing: centered (8), scaled (8)
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 6960, 6959, 6958, 6959, 6959, 6959, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9027401 0.8054708
##
## Tuning parameter 'sigma' was held constant at a value of 0.1
## Tuning
## parameter 'C' was held constant at a value of 10Le realizó la validación cruzada con 6 predictores y 10 pliegues alcanzándose una precisión del 90.43% y un Kappa del 80.86% mostrando una alta concordancia entre las predicciones y los valores reales ajustados por el azar.
11.2.2 N-fold Cross Validation
# Cambia los niveles de la variable de clase 'precio_bin' a nombres válidos
mhp_train_svm$precio_bin <- factor(mhp_train_svm$precio_bin, labels = c("Level0", "Level1"))
# Definimos los controles de entrenamiento para la validación cruzada N-fold
control <- trainControl(method = "cv", number = 5, savePredictions = "final", classProbs = TRUE)
# Creamos el modelo SVM con validación cruzada N-fold
N_F_modelo_svm_cv <- train(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito,
data = mhp_train_svm,
method = "svmRadial",
trControl = control,
preProcess = c("center", "scale"),
tuneGrid = data.frame(sigma = 0.1, C = 10))
# Mostramos los resultados del modelo SVM con validación cruzada N-fold
print(N_F_modelo_svm_cv)## Support Vector Machines with Radial Basis Function Kernel
##
## 7732 samples
## 6 predictor
## 2 classes: 'Level0', 'Level1'
##
## Pre-processing: centered (8), scaled (8)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6185, 6186, 6187, 6185, 6185
## Resampling results:
##
## Accuracy Kappa
## 0.9032597 0.8065112
##
## Tuning parameter 'sigma' was held constant at a value of 0.1
## Tuning
## parameter 'C' was held constant at a value of 10knn_model <- knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")],
k = 5, cl = mhp_train$precio_bin)Le realizó la validación cruzada con 6 predictores y 10 pliegues alcanzándose una precisión del 90.29% y un Kappa del 80.57%.
12 Elección del mejor punto de corte según la curva ROC
Nos disponemos a calcular las curvas ROC para cada uno de los modelo y los visualizamos en un mismo plot para mejor comparabilidad.
# GLM
probs_glm <- predict(glm_mhp, test_glm, type = "response")
# KNN
probs_knn <- attr(prediccion_knn5_test, "prob")
# Decisión Tree
probs_dt <- predict(arbol, mhp_test, type = "prob")[, 2]
# Random Forest
probs_rf <- predict(classifier, mhp_test, type = "prob")[, 2]
# SVM
probs_svm <- predict(modelo_svm, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")], probability = TRUE)
probs_svm <- attr(probs_svm, "probabilities")[, 2]# Calcula las curvas ROC
roc_glm <- roc(mhp_test$precio_bin, probs_glm)
roc_knn <- roc(mhp_test$precio_bin, probs_knn)
roc_dt <- roc(mhp_test$precio_bin, probs_dt)
roc_rf <- roc(mhp_test$precio_bin, probs_rf)
#roc_svm <- roc(mhp_test$precio_bin, probs_svm)
# Configura el gráfico
plot(roc_glm, col = "blue", lwd = 2, legacy.axes = TRUE, main = "Curvas ROC")
lines(roc_knn, col = "red", lwd = 2)
lines(roc_dt, col = "green", lwd = 2)
lines(roc_rf, col = "darkred", lwd = 2)
#lines(roc_svm, col = "orange", lwd = 2)
# Añade leyenda
legend("bottomright", legend = c("GLM", "KNN", "Decisión Tree", "Random Forest", "SVM"),
col = c("blue", "red", "green", "darkred", "orange"), lwd = 2)
CONCLUSIONES
Después de evaluar todos los modelos de aprendizaje automático, se observa que los mejores resultados se obtienen con el modelo de Random Forest y el modelo SVM (Support Vector Machine). Sin embargo, el modelo Random Forest es preferido sobre el modelo SVM por las siguientes razones:
Coste computacional: El modelo Random Forest generalmente tiene un menor costo computacional en comparación con el modelo SVM, especialmente cuando se trabaja con grandes conjuntos de datos. Esto significa que el tiempo de entrenamiento y predicción es más corto y requiere menos recursos computacionales.
Simplicidad del modelo: Aunque ambos modelos son eficaces, el modelo Random Forest es más fácil de entender y de explicar. Los árboles de decisión que forman el Random Forest pueden visualizarse y analizarse fácilmente, lo que facilita la interpretación de los resultados y la identificación de características importantes.
Explicabilidad: La explicabilidad es un aspecto crítico en la ciencia de datos y el aprendizaje automático. Los modelos que son más fáciles de explicar y entender son preferibles en muchas situaciones, especialmente cuando se necesita justificar las decisiones basadas en estos modelos. El modelo Random Forest ofrece una mejor explicabilidad en comparación con el modelo SVM, ya que es más fácil de visualizar y entender cómo se toman las decisiones en el modelo.
Rendimiento similar: A pesar de las ventajas mencionadas anteriormente, la elección entre Random Forest y SVM se basa en gran medida en su rendimiento en la tarea específica. En este caso, ambos modelos proporcionan resultados similares, lo que respalda la decisión de elegir el modelo Random Forest debido a sus otras ventajas.
13 Deep Learning
El Deep Learning es una rama del aprendizaje automático (machine learning) que se enfoca en el entrenamiento de redes neuronales artificiales profundas para aprender y extraer patrones complejos a partir de datos. Estas redes neuronales se componen de múltiples capas de nodos interconectados, lo que les permite aprender representaciones de datos en diferentes niveles de abstracción.
A diferencia de otros métodos de aprendizaje automático, el Deep Learning se caracteriza por su capacidad para aprender automáticamente características relevantes directamente de los datos, sin necesidad de que sean especificadas o extraídas manualmente. Esto lo logra mediante el proceso de entrenamiento, donde la red neuronal ajusta los pesos y los sesgos de sus conexiones para minimizar una función de pérdida y mejorar la capacidad de predicción o clasificación.
El Deep Learning ha mostrado un gran éxito en una amplia variedad de tareas de aprendizaje automático, como reconocimiento de imágenes, procesamiento de lenguaje natural, reconocimiento de voz, traducción automática… Esto se debe a su capacidad para modelar relaciones y patrones complejos en los datos, así como a la disponibilidad de grandes conjuntos de datos y avances en hardware de computación como los GPUs que aceleran el entrenamiento de redes neuronales profundas.
Algunas de las arquitecturas de redes neuronales profundas más utilizadas en Deep Learning incluyen las redes neuronales convolucionales (CNN) para el procesamiento de imágenes, redes neuronales recurrentes (RNN) para el procesamiento de secuencias y las redes neuronales generativas adversariales (GAN) para la generación de contenido nuevo.
Entonces… ¿Cómo funcionan?
Se podría decir que funcionan como un cerebro humano, hay neuronas y cada una de ellas recibe diferentes señales de entrada (𝑥1, 𝑥2, …) de fuentes externas o de otras neuronas. Cada señal de entrada se multiplica por un valor denominado peso que da una idea de la fuerza de dicha conexión (𝑤1, 𝑤2, …). La fase de entrenamiento modifica los pesos, y la neurona calcula una salida, usando Σ y una función de activación 𝑓.
Aplicado a nuestro dataset del mercado inmobiliario en Madrid, vamos a ir entrenando una red neuronal con distintas configuraciones de capas, neuronas y funciones de activación para ver si obtenemos un modelo realmente predictivo del precio de los inmuebles y por consiguiente generar una buena herramienta de clasificación de los mismos.
14 Explicabilidad
15 Series Temporales
Continuando con nuestro proyecto de Machine Learning pero apartándonos bastante del dataset original basado en la predicción del precio y etiquetado de casas en Madrid, procedemos a hacer un estudio de ventas en un famoso centro comercial estadounidense, Walmart.
Para ello hemos obtenido un dataset bastante completo con información de ventas de dicho establecimiento de un repositorio en la siguiente web (kaggle - Walmart Sales Dataset of 45 Stores) con formato long format.
Para el análisis de series temporales conviene saber que son un conjuntos de datos secuenciales que representan observaciones realizadas a lo largo del tiempo, generalmente en intervalos regulares. En este tipo de datos, el orden temporal es fundamental, ya que cada observación está asociada a un momento específico en el tiempo.
El análisis de series temporales se enfoca en comprender y modelar el patrón temporal subyacente en los datos, con el objetivo de predecir o explicar su comportamiento futuro. Las series temporales pueden presentar diferentes características, como tendencias, estacionalidad, ciclos y componentes aleatorios.
Al analizar una serie temporal, se busca identificar patrones y estructuras temporales, como cambios de tendencia, patrones estacionales o efectos de eventos específicos. Esto se logra mediante técnicas estadísticas y modelos matemáticos que permiten modelar y predecir el comportamiento futuro de la serie.
Las series temporales pueden utilizarse en una amplia gama de aplicaciones, como pronósticos económicos, análisis de ventas, previsión del clima, análisis de datos financieros…
En el contexto de nuestro proyecto con el dataset de ventas de Waltmark, nos disponemos a utilizar técnicas de análisis de series temporales para estudiar el comportamiento de las ventas a lo largo del tiempo, identificar patrones estacionales, detectar cambios en las tendencias, analizar el impacto de eventos especiales y realizar pronósticos de las ventas futuras.
15.1 Lectura y preparación de los datos
walmart <- read_csv("walmart.csv")head(walmart, 10) %>%
kbl() %>%
kable_material(c("striped", "hover")) %>%
scroll_box(width = "100%", height = "350px")| Store | Date | Weekly_Sales | Holiday_Flag | Temperature | Fuel_Price | CPI | Unemployment |
|---|---|---|---|---|---|---|---|
| 1 | 05-02-2010 | 1643691 | 0 | 42.31 | 2.572 | 211.0964 | 8.106 |
| 1 | 12-02-2010 | 1641957 | 1 | 38.51 | 2.548 | 211.2422 | 8.106 |
| 1 | 19-02-2010 | 1611968 | 0 | 39.93 | 2.514 | 211.2891 | 8.106 |
| 1 | 26-02-2010 | 1409728 | 0 | 46.63 | 2.561 | 211.3196 | 8.106 |
| 1 | 05-03-2010 | 1554807 | 0 | 46.50 | 2.625 | 211.3501 | 8.106 |
| 1 | 12-03-2010 | 1439542 | 0 | 57.79 | 2.667 | 211.3806 | 8.106 |
| 1 | 19-03-2010 | 1472516 | 0 | 54.58 | 2.720 | 211.2156 | 8.106 |
| 1 | 26-03-2010 | 1404430 | 0 | 51.45 | 2.732 | 211.0180 | 8.106 |
| 1 | 02-04-2010 | 1594968 | 0 | 62.27 | 2.719 | 210.8204 | 7.808 |
| 1 | 09-04-2010 | 1545419 | 0 | 65.86 | 2.770 | 210.6229 | 7.808 |
walmart$Store <- factor(walmart$Store)
walmart$Date <- as.Date(walmart$Date, format = "%d-%m-%Y")
walmart$Holiday_Flag <- as.factor(walmart$Holiday_Flag)summary(walmart)## Store Date Weekly_Sales Holiday_Flag
## 1 : 143 Min. :2010-02-05 Min. : 209986 0:5985
## 2 : 143 1st Qu.:2010-10-08 1st Qu.: 553350 1: 450
## 3 : 143 Median :2011-06-17 Median : 960746
## 4 : 143 Mean :2011-06-17 Mean :1046965
## 5 : 143 3rd Qu.:2012-02-24 3rd Qu.:1420159
## 6 : 143 Max. :2012-10-26 Max. :3818686
## (Other):5577
## Temperature Fuel_Price CPI Unemployment
## Min. : -2.06 Min. :2.472 Min. :126.1 Min. : 3.879
## 1st Qu.: 47.46 1st Qu.:2.933 1st Qu.:131.7 1st Qu.: 6.891
## Median : 62.67 Median :3.445 Median :182.6 Median : 7.874
## Mean : 60.66 Mean :3.359 Mean :171.6 Mean : 7.999
## 3rd Qu.: 74.94 3rd Qu.:3.735 3rd Qu.:212.7 3rd Qu.: 8.622
## Max. :100.14 Max. :4.468 Max. :227.2 Max. :14.313
## str(walmart) ## spc_tbl_ [6,435 × 8] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Store : Factor w/ 45 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Date : Date[1:6435], format: "2010-02-05" "2010-02-12" ...
## $ Weekly_Sales: num [1:6435] 1643691 1641957 1611968 1409728 1554807 ...
## $ Holiday_Flag: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 1 1 ...
## $ Temperature : num [1:6435] 42.3 38.5 39.9 46.6 46.5 ...
## $ Fuel_Price : num [1:6435] 2.57 2.55 2.51 2.56 2.62 ...
## $ CPI : num [1:6435] 211 211 211 211 211 ...
## $ Unemployment: num [1:6435] 8.11 8.11 8.11 8.11 8.11 ...
## - attr(*, "spec")=
## .. cols(
## .. Store = col_double(),
## .. Date = col_character(),
## .. Weekly_Sales = col_double(),
## .. Holiday_Flag = col_double(),
## .. Temperature = col_double(),
## .. Fuel_Price = col_double(),
## .. CPI = col_double(),
## .. Unemployment = col_double()
## .. )
## - attr(*, "problems")=<externalptr>El dataset contiene 6.435 observaciones correspondientes a ventas semanales en cada una de las 45 tiendas Walmart en EEUU entre 2010 y 2012 además de 8 variables (de las cuales son 2 cualitativas y 6 cuantitativas).
A continuación, la descripción de cada una de las variables:
Store: Número de tienda que van del 1 al 45
Date: Día de la semana donde se recoge el dato de ventas. Tiene formato dd-mm-yyyy. Comienza 05-02-2010 hasta 26-10-2012; un total de 143 semanas (2 años y 9 meses)
Weekly_Sales: Ventas totales durante la semana
Holiday_Flag: Si la semana tiene un feriado especial o no. 1 si la semana tiene un feriado 0 si la semana es completamente laboral
Temperature: Temperatura media de la semana donde se producen las ventas
Fuel_Price: Precio del combustible en la región donde se encuentra situada la tienda
CPI: Índice de precios del cliente. Es un indicador utilizado para medir los cambios en el nivel de precios de bienes y servicios a lo largo del tiempo. El CPI es calculado por la Oficina de Estadísticas Laborales (Bureau of Labor Statistics) y es ampliamente utilizado como una medida de la inflación en la economía. Se basa en una “cesta de bienes” que representa los productos y servicios típicos que consume un hogar promedio. Es lo mismo que el IPC (Índice de Precio al Consumo).
Unemployment: Tasa de desempleo en la región donde se encuentra situada la tienda
15.2 Análisis y gráficos univariante
Comenzamos con unas tablas y gráficos simples por variable para ver como se comportan los datos y analizar los estadísticos principales por tienda.
# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}
# Calcular estadísticas de Walmart para la variable Temperature
walmart_stats_temperature <- walmart %>%
group_by(Store) %>%
summarise(min = formatCurrency(round(min(Temperature),2)),
max = formatCurrency(round(max(Temperature),2)),
median = formatCurrency(round(median(Temperature),2)),
mean = formatCurrency(round(mean(Temperature),2)),
sd = formatCurrency(round(sd(Temperature),2)))
# Calcular estadísticas de Walmart para la variable Fuel_Price
walmart_stats_fuel <- walmart %>%
group_by(Store) %>%
summarise(min = formatCurrency(round(min(Fuel_Price),2)),
max = formatCurrency(round(max(Fuel_Price),2)),
median = formatCurrency(round(median(Fuel_Price),2)),
mean = formatCurrency(round(mean(Fuel_Price),2)),
sd = formatCurrency(round(sd(Fuel_Price),2)))
# Calcular estadísticas de Walmart para la variable CPI
walmart_stats_cpi <- walmart %>%
group_by(Store) %>%
summarise(min = formatCurrency(round(min(CPI),2)),
max = formatCurrency(round(max(CPI),2)),
median = formatCurrency(round(median(CPI),2)),
mean = formatCurrency(round(mean(CPI),2)),
sd = formatCurrency(round(sd(CPI),2)))
# Calcular estadísticas de Walmart para la variable Unemployment
walmart_stats_unemployment <- walmart %>%
group_by(Store) %>%
summarise(min = formatCurrency(round(min(Unemployment),2)),
max = formatCurrency(round(max(Unemployment),2)),
median = formatCurrency(round(median(Unemployment),2)),
mean = formatCurrency(round(mean(Unemployment),2)),
sd = formatCurrency(round(sd(Unemployment),2)))
# Mostrar la tabla con los resultados para la variable Temperature
datatable(walmart_stats_temperature, options = list(scrollY = "300px", paging = FALSE), caption = "TEMPERATURE")# Mostrar la tabla con los resultados para la variable Fuel_Price
datatable(walmart_stats_fuel, options = list(scrollY = "300px", paging = FALSE), caption = "FUEL PRICE")# Mostrar la tabla con los resultados para la variable CPI
datatable(walmart_stats_cpi, options = list(scrollY = "300px", paging = FALSE), caption = "CPI")# Mostrar la tabla con los resultados para la variable Unemployment
datatable(walmart_stats_unemployment, options = list(scrollY = "300px", paging = FALSE), caption = "UNEMPLOYMENT")plot1 <- ggplot(walmart, aes(x = "", y = Temperature)) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "ºF", y = "Temperature") +
theme_minimal()
plot2 <- ggplot(walmart, aes(x = "", y = Fuel_Price)) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "$ Gallon", y = "Fuel Price") +
theme_minimal()
plot3 <- ggplot(walmart, aes(x = "", y = CPI)) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "Index", y = "CPI") +
theme_minimal()
plot4 <- ggplot(walmart, aes(x = "", y = Unemployment)) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "Rate", y = "Unemployment") +
theme_minimal()
grid <- grid.arrange(plot1, plot2, plot3, plot4, nrow = 2, ncol = 2)
Nos llama la atención la amplitud de las variables CPI y Unemployment y pasamos a graficarlas por tienda.
plot5 <- ggplot(walmart, aes(x = as.factor(Store), y = CPI, fill = as.factor(Store))) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "Store", y = "CPI", title = "CPI per Store's Region") +
theme_minimal() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
coord_flip()
plot6 <- ggplot(walmart, aes(x = as.factor(Store), y = Unemployment, fill = as.factor(Store))) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "Store", y = "Unemployment", title = "Unemployment Rate per Store's Region") +
theme_minimal() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
coord_flip()
plot5
plot6
Parece que dependiendo de la región donde se encuentre la tienda hay una mayor dispersión del gasto promedio en una cesta de la compra; y la tasa de desempleo también contribuye de manera relevante acentuándo dicha dispersión. El inconveniente es que a través del dataset desconocemos las regiones de EEUU donde se sitúan las 45 tiendas Walmart por lo que no podemos interpretar si tiene o no demasiado sentido, vamos a suponer que sí ya que no es lo mismo una tienda situada en un estado con mayor población que otro, y no es lo mismo una tienda en los estados rurales que en las grandes urbes de ambas costas.
Ahora nos disponemoa a analizar las ventas semanales usando la variable Weekly_Sales.
# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}
# Calcular estadísticas de Walmart
walmart_stats <- walmart %>%
group_by(Store) %>%
summarise(min = formatCurrency(round(min(Weekly_Sales),2)),
max = formatCurrency(round(max(Weekly_Sales),2)),
median = formatCurrency(round(median(Weekly_Sales),2)),
mean = formatCurrency(round(mean(Weekly_Sales),2)),
sd = formatCurrency(round(sd(Weekly_Sales),2)),
sum = formatCurrency(round(sum(Weekly_Sales),2)))
# Mostrar la tabla con los resultados
datatable(walmart_stats, options = list(scrollY = "300px", paging = FALSE), caption = "SALES")# Suma de las ventas semanales por tienda
ventas_por_tienda <- aggregate(Weekly_Sales ~ Store, data = walmart, FUN = sum)
# Orden de las tiendas de menor a mayor ventas
ventas_por_tienda <- ventas_por_tienda[order(ventas_por_tienda$Weekly_Sales, decreasing = FALSE), ]
ventas_por_tienda$Store <- factor(ventas_por_tienda$Store, levels = ventas_por_tienda$Store)
plot7 <- ggplot(ventas_por_tienda, aes(x = Weekly_Sales, y = Store)) +
geom_bar(stat = "identity", fill = "#FFB21C") +
labs(x = "Total Sales (Millions)", y = "Store", title = "Total Sales per Store Aggregated by Sum") +
scale_y_discrete(labels = ventas_por_tienda$Store) +
scale_x_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot7
plot8 <- ggplot(walmart, aes(x = as.factor(Store), y = Weekly_Sales, fill = as.factor(Store))) +
geom_boxplot(fill = "#FFB21C") +
labs(x = "Store", y = "Weekly Sales", title = "Weekly Sales per Store") +
theme_minimal() +
theme(legend.position = "none") +
coord_flip() +
theme(plot.title = element_text(hjust = 0.5)) +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot8
Este boxplot tampoco es muy exclarecedor, en muchas tiendas exite gran dispersión de los datos dentro del rango intercuartílico y hay datos muy atípicos. Esto puede sugerir que la distribución de datos está sesgada. Insistimos en que sin conocer las regiones donde se encuentran las tiendas es muy complicado llegar a más conclusiones o a un mejor entendimiento.
15.3 Análisis y gráficos multivariante
Una vez analizadas todas las variables de forma individual, podemos buscar algunas relaciones entre ellas mediante el siguiente análisis multivariante donde se pueden observar que existen correlaciones muy débiles, sólo pudiendo destacarse la correlación negativa del -0.3 que existe entre el CPI y la tasa de Unemployment; correlación que tiene bastante lógica dado que a mayor tasa de desempleo el gasto promedio en una cesta de la compra se verá redicido. Otro correlación muy débil aunque reseñable es la de 0.18 que existe entre la temperatura y el CPI, quizas a mayor temperatura la población salga más a la calle a consumir… pero es dificilmente demostrable.
variables <- c("Weekly_Sales", "Temperature", "Fuel_Price", "CPI", "Unemployment")
data <- walmart[, variables]
correlation_matrix <- cor(data)
plot9 <- ggcorrplot(correlation_matrix, type = "lower", lab = TRUE, lab_size = 3)
plot9
15.4 Análisis y gráficos de series temporales
Antes de graficar las series temporales de nuestro dataset vamos a sacar por tienda las 3 mayores ventas y el restos de las variables asociadas a las mismas para tratar de identificar si hay algún patrón como que todas se producen en las mismas fechas, que todas tienen unos índeces CPI concretos o que están sujetas a unas tasas de desempleo específicas…
# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}
# Calcular las 5 mayores cifras de ventas para cada tienda
walmart_top5sales <- walmart %>%
group_by(Store) %>%
top_n(3, Weekly_Sales) %>%
summarise(Max_Weekly_Sales = formatCurrency(Weekly_Sales),
Date = Date,
Fuel_Price = formatCurrency(round(Fuel_Price,2)),
CPI = formatCurrency(round(CPI,2)),
Unemployment = formatCurrency(round(Unemployment,2)))
# Mostrar la tabla con los resultados
datatable(walmart_top5sales, options = list(scrollY = "300px", paging = FALSE), caption = "TOP 3 WEEKLY SALES PER STORE")Ahora sacamos las 10 mayores ventas del dataset.
# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}
# Obtener los 10 primeros días con más ventas y aplicar el formato de moneda
top_10_days <- walmart %>%
group_by(Date) %>%
summarise(total_sales = sum(Weekly_Sales)) %>%
top_n(10, total_sales) %>%
arrange(desc(total_sales)) %>%
mutate(total_sales = formatCurrency(total_sales))
datatable(top_10_days, options = list(scrollY = "300px", paging = FALSE), caption = "TOP 10 DAYS WITH HIGHEST SALES")Se observa claramente como en ambas tablas las semanas de Navidad y Acción de Gracias en noviembre tienen el record de ventas todos los años.
GRÁFICO SERIE TEMPORAL
Comenzamos realizando con unos gráficos de líneas temporales en las variables más relevantes que deseamos analizar.
walmart <- arrange(walmart, Date)
plot10 <- ggplot(walmart, aes(x = Date, y = Weekly_Sales)) +
geom_line() +
labs(x = "Date", y = "Weekly Sales", title = "Weekly Sales Time Series Plot") +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6)) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5))
plot10
plot11 <- ggplot(walmart, aes(x = Date, y = Temperature)) +
geom_line() +
labs(x = "Date", y = "Temperature", title = "Temperature Time Series Plot") +
theme(plot.title = element_text(hjust = 0.5)) +
theme_minimal()
plot12 <- ggplot(walmart, aes(x = Date, y = Fuel_Price)) +
geom_line() +
labs(x = "Date", y = "Fuel Price", title = "Fuel Price Time Series Plot") +
theme(plot.title = element_text(hjust = 0.5)) +
theme_minimal()
plot13 <- ggplot(walmart, aes(x = Date, y = CPI)) +
geom_line() +
labs(x = "Date", y = "CPI", title = "CPI Time Series Plot") +
theme(plot.title = element_text(hjust = 0.5)) +
theme_minimal()
plot14 <- ggplot(walmart, aes(x = Date, y = Unemployment)) +
geom_line() +
labs(x = "Date", y = "Unemployment", title = "Unemployment Time Series Plot") +
theme(plot.title = element_text(hjust = 0.5)) +
theme_minimal()
grid <- grid.arrange(plot11, plot12, plot13, plot14, nrow = 2, ncol = 2)
El gráfico de la temperatura muestra un comportamiento totalmente normal donde se pueden identificar a simple vista las estaciones y sus variaciones de temperatura entre invierno y verano.
El gráfico del precio del combustible muestra una relativa estabilidad en precios sufriendo una subida a partir del primer trimestre del 2010 teniendo su punto más álgido en el segundo trimestre del año; punto a partir del cual vuelve a bajar situándose cerca de los 2,5 dólares por galón en el tercer trimestre del año. No es así en el cuarto trimestre donde se origina una vertiginosa subida del precio por galón hasta el segundo trimestre del 2011, pasando de 2,5 a 4 dólares. Desde este punto hasta el final del año se produce otra severa caída de precios cerrando el año 2011 alrededor de los 3,25 dólares. Durante el año 2012 se van produciendo unas subidas y bajadas de precios con bastante apuntamiento entre los 3,25 y 4 dólares.
El gráfico del CPI muestra un comportamiento monótono creciente pero de manera apaisada que nos indica que el gasto medio en la lista de la compra va aumentando lentamente de manera progresiva. No obstante el gráfico muestra gran varianza en este índice debido precisamente a lo ya comentado con anteriorida; hay tiendas Walmart situadas en regiones de EEUU donde índice de población y el poder adquisitivo es menor y otras donde es mayor, diferencias entre estados rurales y grandes urbes.
El gráfico de la tasa Unemployment parece que se muestra estable a lo largo del 2010 y 2011 hasta comenzar 2012 donde generalmente se va produciendo una caída en las tasas de desempleo.
15.4.1 Estacionariedad
Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada. Dadas las características de la serie realizamos una descomposición STL para corroborar que existe tendencia y estacionalidad.
La descomposición STL (Seasonal and Trend decomposition using Loess) es un método utilizado para descomponer una serie temporal en sus componentes de tendencia, estacionalidad y residuos. Utiliza un enfoque no paramétrico basado en el ajuste de suavizadores locales (Loess) para estimar la tendencia y la estacionalidad de la serie temporal. La tendencia representa el patrón general y de largo plazo en los datos, mientras que la estacionalidad captura los patrones recurrentes y periódicos. Los residuos son la parte de la serie temporal que no puede ser explicada por la tendencia y la estacionalidad.
# Convertir el dataset a un tsibble
walmart_tsibble <- as_tsibble(walmart, key = Store, index = Date)walmart_tsibble_train <- walmart_tsibble %>% filter_index(. ~ "2012-06-29")
walmart_tsibble_test <- walmart_tsibble %>% filter_index("2012-07-06" ~ .)walmart_tsibble_train %>%
model(seats = feasts:::STL(Weekly_Sales)) %>%
components() %>%
autoplot() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6)) +
theme(panel.grid.major = element_line(color = "grey"),
panel.grid.minor = element_line(color = "grey"),
panel.border = element_rect(color = "grey", fill = NA))
Se confirman visualmente las estacionalidades supuestas.
walmart_tsibble_train %>%
model(seats = feasts:::STL(Fuel_Price)) %>%
components() %>%
autoplot() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_line(color = "grey"),
panel.grid.minor = element_line(color = "grey"),
panel.border = element_rect(color = "grey", fill = NA))
walmart_tsibble_train %>%
model(seats = feasts:::STL(CPI)) %>%
components() %>%
autoplot() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_line(color = "grey"),
panel.grid.minor = element_line(color = "grey"),
panel.border = element_rect(color = "grey", fill = NA))
walmart_tsibble_train %>%
model(seats = feasts:::STL(Unemployment)) %>%
components() %>%
autoplot() +
theme(legend.position = "none") +
theme(plot.title = element_text(hjust = 0.5)) +
theme(panel.grid.major = element_line(color = "grey"),
panel.grid.minor = element_line(color = "grey"),
panel.border = element_rect(color = "grey", fill = NA))
Estacionariedad en MEDIA
walmart_tsibble %>% features(Weekly_Sales, unitroot_kpss)## # A tibble: 45 × 3
## Store kpss_stat kpss_pvalue
## <fct> <dbl> <dbl>
## 1 1 0.476 0.0471
## 2 2 0.0787 0.1
## 3 3 0.608 0.0219
## 4 4 0.920 0.01
## 5 5 0.672 0.0161
## 6 6 0.0365 0.1
## 7 7 0.512 0.0390
## 8 8 0.187 0.1
## 9 9 0.669 0.0163
## 10 10 0.131 0.1
## # … with 35 more rowsEstacionariedad en VARIANZA
15.4.2 Determinación del modelo
15.4.2.1 Modelo A
15.4.2.2 Modelo B
15.4.3 Fase de estimación y contraste
15.4.4 Diagnosis de residuos
Test de Media = 0
Test de Autocorrelaciones
Test de Homocedasticidad
Test de Normalidad
15.4.5 Fase de predicción
Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones y calculamos sus errores y los comparamos con los residuos.